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Introduction 


It  is  clear  from  both  a  logistical  and  an  economic  viewpoint  that  the  combustion  of 
hydrocarbon  fuels  will  continue  to  play  a  central  role  in  US  Air  Force  operations.  In 
addition,  with  world  events  demanding  enhanced  flexibility  in  the  sources  of  fuels,  it  is 
becoming  increasingly  important  to  improve  scientific  knowledge  of  the  combustion 
properties  of  different  fuel  blends.  Chemically  reacting  flows  utilizing  hydrocarbon  fuel 
blends  occur  in  a  variety  of  energy  conversion  processes  such  as  combustion,  propulsion, 
and  fuel  reforming,  to  name  just  a  few.  They  are  also  relevant  to  many  material  synthesis 
technologies.  If  the  chemical  transformations  and,  in  some  cases,  attendant  energy  release 
can  be  made  to  happen  in  a  well-controlled  and  well-defined  fashion,  the  goal  of  either 
liberating  heat,  partially  oxidizing  a  fuel,  generating  electrical  power,  or  synthesizing 
advanced  materials  should  be  achievable  with  high  efficiency  and  minimal  pollution. 

Hydrocarbon  fuels  will  remain  a  major  source  of  energy  well  into  the  second  half  of  the 
21st  century  and,  despite  dire  warnings  about  their  limited  supply,  known  resources  have 
actually  increased  over  the  past  decade.  Nevertheless,  finite  supplies  will  continue  to 
exert  pressure  on  the  efficient  use  of  these  fuels,  especially  as  the  price  of  oil  continues  to 
skyrocket.  In  the  engineering  of  chemically  reacting  flows,  increased  efficiency  and 
reduced  pollution  can  be  achieved  via  an  integrated  approach  that  extends  the  state-of- 
the-art  of  both  experimental  and  computational  methodologies.  By  making  advances  in 
each  of  these  areas  and  by  integrating  them  in  well-conceived  research  programs, 
scientists  will  be  able  to  have  a  dramatic  impact  on  the  design  of  technologies  involving 
energy  conversion  and  combustion.  One  of  the  most  technically  challenging  engineered 
systems  of  importance  to  the  Air  Force  in  which  chemically  reacting  flows  play  a  critical 
role  are  gas  turbines  (GTs).  For  aeropropulsion  applications,  there  are  no  alternative 
energy  replacements  of  GTs  in  sight.  In  addition,  the  majority  of  the  electrical  power  to 
be  added  in  the  United  States  and  around  the  world  through  2015  will  be  based  on  GTs. 
Advances  in  GT  engineering  would  inevitably  affect  the  entire  combustion  industry,  and 
the  economic  repercussions  of  these  advances  would  be  amplified  further  as  critical 
combustion  issues  for  GTs  are  also  important  to  the  transportation  industry.  Moreover, 
the  economic  payoff  for  US  Military  Operations  could  be  enormous. 

Many  contributions  are  needed  to  advance  the  frontiers  of  the  science  behind  such 
engineered  systems,  and  thereby  to  enhance  the  nation's  economic  base  and  help  stabilize 
it  against  foreign  competition  and  dependencies.  The  research  will  focus  on  the 
development  of  advanced  computational  methodologies  that  will  enable  reacting  flow 
simulations  which  are  more  rapid  and  more  accurate  than  those  currently  feasible,  and 
which,  on  both  counts,  are  capable  of  deepening  an  understanding  of  the  fluid  dynamics 
and  aerothermochemistry  underlying  many  vital  technologies.  Specifically,  this  research 
has  considered  numerical  algorithms  designed  for  the  solution  of  gas-phase  combustion 
with  detailed  transport  and  finite-rate  chemistry.  To  help  achieve  these  goals,  a 
companion  experimental  program  has  been  initiated  in  which  the  complexity  of  the 
various  systems  is  being  dissected  into  well-defined  laboratory-scale  problems,  from 
which  data  can  be  provided  for  the  validation  of  the  computational  models. 


Overview  of  the  Implicit-Compact  Solver 

The  implicit-compact  methods  studied  in  this  granting  period  have  been  designed  to  meet 
two  well-known  challenges  in  modeling  time-dependent  combustion:  the  stiffness 
induced  by  the  vastly  disparate  timescales  in  the  chemistry,  which  calls  for  implicit  time 
integration;  and  the  significant  spatial  structure  in  the  flow  field,  which  cannot  be 
captured  without  high  resolution,  low  diffusivity  spatial  discretizations.  The  main  idea  of 
a  compact  scheme  discretization  is  to  construct  algebraic  relationships  between  the  values 
of  a  function  and  of  its  derivative  at  the  nodes  of  a  grid.  These  equations  are  written  in 
matrix  form;  the  matrices  are  banded,  generally  tridiagonal,  and  hence  can  be  inverted 
efficiently.  The  coefficients  are  constant  for  a  given  grid  and  are  defined  by  matching  the 
terms  of  Taylor  series  expansions.  The  spatial  discretizations  used  in  this  work  have  a 
variable  order  of  accuracy,  depending  on  the  grid  spacing  and  the  presence  of  steep 
gradients  near  the  domain  boundaries.  The  maximum  order  is  six  and  the  minimum  is 
three.  Quite  apart  from  their  classical  order  of  accuracy,  the  particular  advantage  of  the 
compact  schemes  is  their  “spectral-like”  resolution  of  moderately  high  wave  numbers.  In 
practice,  this  allows  good  accuracy  over  a  long  time  with  many  fewer  grid  points  than 
would  be  required  by  a  traditional  low  order  finite  difference  method. 

In  the  implicit-compact  solver,  the  governing  partial  differential  equations  (PDEs)  are 
semi-discretized  using  a  compact  finite  difference  procedure  (a  finite  volume  method 
could  be  used  too).  Then,  after  the  spatial  discretization,  the  system  of  ordinary 
differential  equations  is  discretized  with  an  A-stable  backward  difference  formula  (BDF), 
following  the  “method  of  lines”  approach.  The  resulting  nonlinear  algebraic  system  is 
solved  by  a  damped  inexact  Newton’s  method.  An  approximate  solution  to  the  linearized 
Newton  system  is  obtained  using  an  iterative  Krylov  method  (GMRES)  with  an 
appropriate  preconditioner  (incomplete  LU  decomposition  with  a  scaling/reordering 
preprocessor).  The  solver  has  been  thoroughly  tested  on  problems  with  known  analytical 
solutions,  thus  verifying  the  correctness  of  all  temporal  and  spatial  discretizations.  Seeing 
as  the  preconditioner  is  the  most  expensive  part  of  the  solution  process,  coarse-grain 
parallelism  was  introduced  into  this  module  by  means  of  restricted  additive  Schwarz 
domain  decomposition,  implemented  in  a  shared  memory  context  by  OpenMP  pragmas. 
The  payoff  here  was  not  significant  (due  to  memory  bandwidth  issues,  the  code  could 
only  run  effectively  on  4-6  threads  at  once),  and  the  extension  of  the  domain 
decomposition  algorithm  to  a  distributed  memory  context  remains  a  task  for  future  work. 
As  a  result  of  an  AFOSR  DURIP  award,  the  hardware  for  this  work  is  already  in  place:  a 
cluster  with  128  cores,  512  GB  of  RAM,  and  several  terabytes  of  disk,  all  connected  via  a 
high-speed  DDR  Infiniband  fabric. 


The  code  has  been  written  in  Fortran.  By  now,  almost  all  of  the  modules  have  been  (at 
least  partly)  modernized  to  take  advantage  of  Fortran  90/95  enhancements  in  the  area  of 
array  processing,  user-defined  data  structures,  and  data  encapsulation  in  modules.  The 
only  exception  is  the  MC64  code,  which  has  been  taken  from  the  Flarwell  Subroutine 
Library.  The  Fortran  for  all  of  the  problem-specific  subroutines  has  been  produced  by  a 
Mathematica  code  generation  tool,  as  noted  below.  Partial  listings  of  this  script  and  of  the 
residual  and  Jacobian  subroutines  created  with  it  are  attached  as  Appendices. 


Accomplishments:  Modeling 


(a)  Computing  the  pressure  field :  Originally,  the  plan  was  to  work  with  a  velocity- 
vorticity  formulation  of  the  fluid  dynamical  problem,  a  decision  based  primarily  on  the 
experience  obtained  in  using  vorticity-based  methods  to  model  flames.  A  side  benefit 
here  was  that  there  would  be  no  need  to  contend  with  the  signature  challenge  of  low 
speed  flows:  computing  the  pressure  field.  Before  long,  however,  it  became  clear  that  the 
velocity-vorticity  formulation  leads  to  intractable  linear  systems  in  the  Newton  iteration 
for  time-dependent  problems  discretized  in  space  with  compact  schemes.  Accordingly, 
during  the  first  six  months  of  the  grant,  the  intended  approach  to  the  fluid  dynamics  was 
abandoned  in  favor  of  a  primitive  variables  formulation.  Since  acoustic  effects  were  not 
of  primary  interest,  the  zero  Mach  number  approximation  was  employed  in  formulating 
the  governing  equations.  As  is  standard  in  the  modeling  of  low-speed  flows,  the 
hydrodynamic  pressure  was  taken  as  the  independent  variable,  thermodynamic  pressure 
was  constant,  and  the  density  was  recovered  from  the  ideal  gas  law.  Note  that,  with  a 
fully  coupled  solver,  continuity  can  be  enforced  even  though  dynamic  pressure  is  retained 
as  the  corresponding  unknown.  Also,  by  constructing  the  numerical  method  to  take 
advantage  of  the  natural  coupling  of  the  variables,  a  Poisson  equation  for  pressure  and  a 
pressure  projection  step  are  not  needed.  In  studying  flow  in  a  pipe,  it  was  initially 
observed  that  the  smoothness  of  the  solutions  tended  to  be  very  sensitive  to  the  choice  of 
grid;  specifically,  it  was  found  that  local  under-resolution  of  the  radial  velocity  field  led 
to  noticeable  oscillations  in  the  pressure  field  along  most  of  the  length  of  the  pipe.  The 
lesson  learned  was  that  the  tight,  global  coupling  of  all  unknowns  in  the  implicit-compact 
discretization  places  high  demands  on  the  adequate  spatial  resolution  of  all  structures  and 
boundary  layers  in  the  flow  field.  Once  this  was  achieved,  the  solver  delivered  strikingly 
smooth,  accurate  solutions,  the  likes  of  which  it  was  impossible  to  achieve  with 
traditional  low  order  methods.  These  results  have  demonstrated  that  the  notorious 
“pressure-velocity  decoupling”  problem,  which  manifests  itself  as  an  oscillatory  pressure 
field  in  finite  difference  calculations  using  collocated  grids,  can  be  overcome  by  using 
compact  schemes  to  discretize  in  space: 
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Contours  of  dynamic  pressure,  compact  scheme  (color)  vs.  low  order  method  (black  lines). 
Left:  steady  pipe  flow  (Re=500,  pipe  radius=0.4cm).  Right:  oscillating  cold  jet  flow  at  0.025  s. 
The  small-amplitude  grid-scale  noise  in  the  low  order  solutions  is  due  to  the  pressure-velocity 
decoupling. 


(b)  Unsteady  multicomponent  flows  with  thermal  mixing :  Before  attempting  to  compute 
reacting  flows  it  was  imperative  to  move  beyond  simple  flows  and  demonstrate  the 
capability  of  the  implicit-compact  solver  to  model  complicated  multicomponent  flows 
with  thermal  and  species  mixing.  The  work  focused  on  jet  flows  in  quasi-open 
axisymmetric  geometries,  where  the  flow  configuration  was  chosen  to  mimic  that  of  the 
diffusion  flames  which  are  the  goal  of  this  research.  Here,  “quasi-open”  means  that  a 
solid  wall  was  placed  in  the  radial  far-field  of  the  computational  domain,  at  a  distance  of 
many  jet  radii  from  the  centerline.  Though  this  introduced  a  flow  recirculation  zone  near 
this  wall,  the  velocities  were  very  small  and  any  difficulty  of  computing  the  recirculation 
was  more  than  compensated  by  the  ability  to  use  a  simple  “no  slip”  Dirichlet  boundary 
condition  there  for  the  velocity  field.  This  took  the  complicated  issue  of  open  boundary 
conditions  off  the  table  where  it  was  most  likely  to  be  a  sticking  point,  in  the  radial  far- 
field  (homogeneous  Neumann  outflow  conditions  are  probably  feasible  for  the  high  order 
solver  since  the  grids  used  in  flame  calculations  are  typically  very  long  in  the  axial 
direction).  Calculations  of  an  oscillating  heated  jet  issuing  into  quiescent  cool  air  were 
undertaken.  The  reference  calculations  used  the  same  low  order  discretization  employed 
in  a  previously  developed  flame  code,  namely,  second  order  centered  differences  for 
diffusive  terms  and  first  order  (monotone)  upwinding  for  convective  terms.  These 
solutions  were  ruined  by  strong  artificial  viscosity  from  the  upwind  differencing.  By 
contrast,  the  calculations  with  the  compact  scheme  discretization  successfully  captured 
the  spatial  structure  of  the  flow,  revealing  a  marked  “pinching”  in  the  temperature  field 
that  was  qualitatively  similar  to  the  thermal  structure  of  the  time-varying  diffusion  flame 
studied  at  Yale  in  the  laboratory  of  Marshall  Long.  Both  calculations  were  run  with  a 
time  step  small  enough  to  ensure  that  the  spatial  error  dominated.  A  comparison  of  the 
results  is  presented  in  the  figure  below. 


Temperature  contours  in  a  forced,  heated  jet  flow  (Re=500,  AT=100K),  compact  scheme  (CS) 
vs.  low  order  (LO)  solution.  The  forcing  frequency  is  20  Hz.  The  images  are  taken  at  0.0, 
0.0125,  0.025,  and  0.0375  s. 
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(c)  Complex  chemistry  capability :  A  model  premixed  flame  was  studied  in  order  to 
isolate  the  numerical  challenges  of  realistic  combustion  thermochemistry  from  those 
relating  to  the  fluid  dynamics  (the  flow  field  is  imposed  in  the  model).  Three  convection- 
diffusion-reaction  equations  with  exponentially  nonlinear  two-step  Arrhenius  chemistry 
were  solved  for  temperature  and  two  reacting  species  in  a  two-dimensional  axisymmetric 
geometry.  These  calculations  have  provided  a  proof-of-concept  that  the  implicit-compact 
methods  can  deliver  accurate  and  efficient  solutions  to  stiff  multi-step  chemistry 
combustion  problems.  These  types  of  problems  arise  in  modeling  the  combustion  of 
aviation  fuels  of  interest  to  the  Air  Force,  and  they  create  significant  numerical 
challenges  for  both  explicit  methods  (CFL  restrictions)  and  many  splitting  methods 
(additional  accuracy  limitations  due  to  splitting  errors).  Calculations  of  the  steady,  two- 
dimensional  model  flame  were  performed  with  both  the  compact  scheme  and  the  low 
order  semi-discretizations,  and  the  steady  solutions  were  then  used  to  initialize  time- 
dependent  calculations,  where  the  imposed  flow  was  periodically  varying  (much  like  in 
the  oscillating  jet  flow  discussed  above).  In  the  transient  simulation,  comparison  of  the 
compact  scheme  solution  with  the  low  order  solution  revealed  the  presence  of  numerical 
dif  fusion  in  the  latter.  This  can  be  observed  most  clearly  in  the  damping  of  the  temporal 
oscillation  of  the  intermediate  species.  This  artifact  of  the  spatial  discretization  had 
hampered  earlier  low  order  simulations  of  time-varying  diffusion  flames.  The  implicit- 
compact  solver  succeeds  in  generating  spatially  accurate  solutions  for  such  problems,  and 
since  it  allows  for  large  time  steps  the  computational  cost  is  reasonable. 
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Evolution  of  intermediate  species  concentration  in  the  model  premixed  flame  problem, 
compact  scheme  (CS)  vs.  low  order  (LO)  solution.  The  frequency  of  oscillation  is  10  Hz.  The 
images  are  taken  at  0,0,  0.02,  0.04,  and  0.06  s,  respectively. 

(d)  Experimental  validation  of  the  solver.  Concurrent  with  the  numerical  development 
was  the  construction  of  the  experimental  configuration  for  validating  the  new  implicit- 
compact  solver,  which  will  be  essential  as  the  work  transitions  to  detailed  chemistry 
flame  simulations.  Toward  the  end  of  the  granting  period,  a  burner  was  tested  in  which 


fuel  flows  from  a  0.4  cm  inner  diameter  vertical  tube  (wall  thickness  0.038  cm)  into  a 
concentric,  7.4  cm  diameter  oxidizer  coflow. 


CAD  drawing  of  the  forced-flow  burner  used  in  the  time-varying  diffusion  flame  experiments. 


A  speaker  in  the  plenum  of  the  fuel  jet  allows  a  periodic  perturbation  to  be  imposed  on 
the  exit  parabolic  velocity  profile.  The  fuel  is  diluted  with  nitrogen  and  the  velocities 
have  been  carefully  tuned  so  as  to  produce  a  flame  with  negligible  soot  that  is  lifted 
above  the  burner  surface,  preventing  heat  transfer  from  the  flame  to  the  burner.  These 
operating  conditions  have  been  designed  to  simplify  the  calculations  in  such  a  way  that 
physical  realism  is  not  sacrificed.  Recently,  laser  diagnostic  techniques  such  as  Raman 
and  Rayleigh  scattering  have  been  used  to  measure  temperature  and  major  species 
profiles  from  flames  using  this  burner. 


Accomplishments:  Numerical  Methods  and  Software 

(a)  Efficient  Jacobian  operations  for  compact  scheme  discretizations :  The  first  hurdle  in 
developing  the  implicit-compact  methodology  has  been  to  devise  a  set  of  highly  efficient 
numerical  algorithms  for  the  formation,  multiplication,  and  element-extraction  of  the 
Jacobian  matrix,  and  to  develop  the  means  to  implement  these  algorithms  effectively  in  a 
large  scientific  code.  The  methods  used  to  accomplish  these  tasks  have  come  to  maturity 
over  the  past  few  years.  Though  dense  by  any  typical  criterion  of  matrix  sparsity, 
compact  scheme  Jacobians  possess  a  latent  structure  based  on  the  fact  that  they  arise  from 
a  discretization  of  a  PDE.  The  idea  of  decomposing  the  Jacobian  into  “local”  and 
“spatial”  components  exploits  this  structure  to  achieve  a  form  of  data  compression.  To 
illustrate  the  idea:  if  U  is  the  independent  variable  and  F  =  F(U,Ux)  is  the  residual,  which 
is  a  function  of  U  and  its  spatial  derivative  Ux,  the  Jacobian  matrix  in  Newton’s  method 
can  be  reconstructed  by  considering  Ux  to  be  independent  of  U,  forming  the  “local” 
Jacobian  components  dF/dU  and  the  “spatial”  Jacobian  components  dF/d(Ux),  and 
observing  that  J  =  dF/dU  +  (dF/d(Ux))*Dx,  where  Dx  is  the  appropriate  coefficient  of  the 
differentiation  matrix  used  to  calculate  Ux  from  U.  The  result  of  using  this  decomposition 


is  that,  even  for  the  most  memory  intensive  combustion  problems,  a  Jacobian  arising 
from  a  compact  scheme  spatial  discretization  can  be  stored  with  less  memory  than  a 
conventional  finite  difference  method  would  need  to  store  an  equivalent  low  order 
Jacobian.  Moreover,  once  generated,  the  compact  scheme  Jacobian  can  be  applied  to  a 
vector  extremely  efficiently,  since  J*V  =  (dF/dU)*V  +  (dF/d(Ux))*(Dx*V)  =  (dF/dU)*V 
+  (dF/d(Ux))*  Vx.  Here  it  is  seen  that  J*V  is  equal  to  the  sum  of  two  dot  products,  one  of 
the  “local  components”  with  the  vector  V  and  the  other  of  the  “spatial  components”  with 
the  vector  Vx.  The  fundamental  operation  in  any  iterative  linear  algebra  routine  is  the 
matrix-vector  product.  With  this  algorithm  it  can  be  computed  in  O(N)  flops. 

(b)  Linear  algebra  enhancements  /.  Robust  iterative  methods :  A  fully  coupled  solution 
paradigm  makes  little  sense  without  robust  and  efficient  numerical  linear  algebra 
algorithms.  These  already  exist  for  the  low  order  spatial  discretizations  commonly  used 
in  computational  combustion;  for  the  compact  discretizations,  they  have  been  discovered 
and  fine-tuned  during  the  tenure  of  this  grant.  In  the  first  year,  in  the  course  of  studying 
pipe  flows,  it  was  found  that  the  use  of  large  grids  (i.e.,  many  points)  led  to  stagnation  in 
the  preconditioned  Bi-CGSTAB  linear  solver.  The  difficulties  were  more  severe  than 
anything  encountered  using  Bi-CGSTAB  in  solving  combustion  problems  within  the  past 
fifteen  years.  To  overcome  them,  a  new  iterative  linear  algebra  module  was  developed 
and  integrated  into  the  implicit-compact  solver.  Based  on  GMRES,  it  enjoys  the 
monotonicity  and  enhanced  robustness  of  this  method,  while  retaining  the  state-of-the-art 
MC64-ILUT  preconditioner  that  has  proven  effective  for  the  challenging  linear  systems 
produced  by  compact  semi-discretization.  The  implementation  made  some  sacrifices  of 
efficiency  for  greater  reliability,  such  as  eschewing  restarts  and  performing  modified 
Gram-Schmidt  with  full  re-orthogonalization  for  the  Amoldi  process,  though  the  penalty 
incurred  was  minimal  for  time-dependent  problems,  where  the  Jacobian  matrices  have 
large  diagonal  terms  and  the  linear  system  solution  is  fast.  In  any  case,  even  in  the 
absence  of  optimizations,  the  new  linear  solver  not  only  made  it  possible  to  solve 
problems  that  had  previously  caused  the  code  to  fail  (e.g.,  a  large  binary  mixing  problem 
on  very  nonuniform  grids),  but  also  displayed  significant  performance  gains  over  its 
predecessor  (e.g.,  showing  between  a  two-  and  five-fold  decrease  in  the  time  spent 
solving  linear  systems).  A  decade  ago,  memory  constraints  made  the  use  of  GMRES  in 
flame  calculations  much  more  difficult  than  it  is  today. 

(c)  Linear  algebra  enhancements  II.  Robust  preconditioning :  The  basic  approach  of 
applying  a  purely  algebraic  MC64-ILUT  preconditioning  algorithm  to  a  pre-sparsified  or 
“partial”  Jacobian  matrix  is  sound,  but  also  potentially  expensive  and  difficult  to  optimize 
due  to  the  parameters  in  the  algorithm.  Early  on,  by  performing  a  sequence  of  tests  on 
basic  fluid  dynamics  test  problems  one  could  sample  enough  of  the  parameter  space  of 
the  preconditioner  to  come  to  a  satisfactory  understanding  of  how  at  least  to  get  the 
incomplete  factorization  to  work.  With  ILUT,  it  became  clear  that  the  drop  tolerance  was 
the  more  expensive  and  less  effective  parameter  to  tune.  Hence,  the  strategy  originally 
employed  was  to  compute  with  a  modest  fill-in  parameter  and  to  recompute  with  more 
fill  if  the  linear  solver  failed.  As  long  as  this  worked,  it  was  fine;  however,  as  the 
difficulty  of  the  physical  models  and  the  size  of  the  problems  were  scaled  up,  both  the 
heightened  importance  of  the  preliminary  sparsification  phase  and  the  pressure,  for  the 
sake  of  computational  efficiency,  to  form  the  preconditioner  less  frequently  clarified  a 


new  question:  when  the  linear  solver  fails,  is  this  due  to  a  specific  problem  with  the  1LU 
or  rather  simply  to  an  out-of-date  or  otherwise  ineffective  preconditioner?  In  the  effort  to 
answer  this  question,  the  linear  solver  module  was  enhanced  by  the  addition  of  a  number 
of  inexpensive  sanity  checks  that  helped  us  to  interpret  the  results  of  successful  linear 
solves  and  to  diagnose  the  cause  of  failures.  By  far  the  most  frequent  cause  of  failure  was 
instability  in  the  ILU.  This  was  deduced  by  comparing  the  condition  estimates  of  the  ILU 
and  of  the  preconditioned  linear  system,  or  rather,  of  the  Hessenberg  matrix  constructed 
by  GMRES,  which  represents  the  projection  of  this  linear  system  onto  the 
orthonormalized  Krylov  basis.  A  poorly  conditioned  ILU  indicates  instability  in  the 
incomplete  factorization  process;  given  a  well  conditioned  ILU,  a  poorly  conditioned 
Hessenberg  suggests  that  the  preconditioner  is  inaccurate.  The  diagnosis  is 
straightforward,  and  in  the  former  case,  so  is  the  cure:  perturb  the  diagonal  of  the  partial 
Jacobian  and  recomputed  the  ILU.  This  understanding  was  a  breakthrough  for  two 
reasons.  First,  it  meant  that  many,  perhaps  the  majority  of,  failed  linear  solves  could  be 
rectified  at  essentially  no  increase  in  computational  cost.  Applying  a  diagonal 
perturbation  is  a  negligible  expense  compared  to  the  ILU,  the  ILU  can  be  performed 
without  increasing  the  level  of  fill-in,  and  if  the  perturbation  is  not  large  the  resulting 
preconditioned  Krylov  process  can  still  converge  in  a  reasonable  number  of  iterations.  By 
contrast,  the  previous  approach  to  addressing  failure  in  the  linear  solver  (recomputing  the 
ILU  with  more  fill-in)  has  a  cost  in  terms  of  both  memory  usage  and  time  which  grows 
unpredictably  with  increasing  fill.  Second,  the  use  of  diagonal  perturbations  has  allowed 
stabilized  ILU  factorizations  of  compact  scheme  partial  Jacobians  at  much  larger  time 
steps  than  otherwise  possible.  This  has  again  made  it  feasible  to  calculate  steady  solutions 
with  the  compact  scheme  solver  for  some  large-scale  problems  where  previously  this  had 
been  out  of  the  question  (see  the  figure  below,  and  discussion  in  (d)). 


(d)  Constructing  good  initial  conditions  for  time-dependent  problems:  The  semi- 
discretized  governing  equations  for  a  transient  flame  problem  in  the  primitive  variable 
formulation  are  a  system  of  differential-algebraic  equations  (DAEs).  As  is  well  known,  a 
DAE  requires  consistent  initialization:  it  is  important  to  construct  an  initial  condition  that 
satisfies  the  algebraic  constraints  so  as  to  avoid  boundary  layers  and  maintain  accuracy  in 
the  early  phase  of  the  computations.  Moreover,  if  high  order  spatial  discretizations  were 
used  in  moving  from  the  original  PDE  to  the  DAE,  this  initial  condition  must  itself  meet 
strict  smoothness  requirements.  A  basic  challenge  for  those  employing  high  order 


methods  is  how  to  generate  such  an  initial  condition.  In  practice,  one  often  wishes  to  use 
a  steady-state  solution  for  this  purpose.  This  was  the  approach  assumed  at  the  inception 
of  this  research  program.  However,  it  was  very  difficult  to  generate  a  steady  compact 
scheme  solution  from  scratch  using  Newton  iteration  with  pseudo-transient  continuation, 
the  method  of  choice  for  problems  discretized  with  traditional  low  order  finite 
differences.  The  source  of  the  difficulty  is  that,  contrary  to  their  name  and  their 
reputation,  the  effective  computational  stencils  of  compact  schemes  are  very  wide  (the 
inverse  of  a  tridiagonal  matrix  is  structurally  dense),  so  Jacobian  matrices  based  on 
compact  scheme  discretizations  can  be,  and  at  steady  state  generally  are,  very  far  from 
diagonal  dominance,  and  hence  very  difficult  to  precondition  using  known  methods.  As  a 
workaround,  a  low-cost  and  effective  means  of  computing  satisfactorily  smooth  and 
accurate  initial  conditions  for  high  order  flow  simulations  was  developed.  The  method 
takes  the  '‘inadequate”  steady-state  solutions  produced  by  low  order  solvers  and  marches 
them  in  pseudo-time  with  moderately  large  time  steps  using  the  implicit-compact  solver. 
This  procedure  has  been  greatly  facilitated  by  progress  in  learning  how  to  stabilize  the 
1LU  preconditioner  at  such  time  steps,  as  described  above. 

(e)  Code  generator :  At  present,  the  entire  computational  kernel  of  the  implicit-compact 
solver  (the  residual  function,  the  Jacobian  operations,  and  all  of  the  problem-specific 
preconditioner  routines)  can  be  generated  automatically  in  Fortran  or  C,  by  simply 
entering  the  partial  differential  equations  and  the  boundary  conditions  in  standard 
mathematical  notation  and  executing  a  Mathematica  script.  For  simplicity,  vector 
calculus  operators  (Grad,  Div)  can  be  used  in  formulating  the  problem,  which  another 
freely  available  Mathematica  package  translates  into  the  correct  partial  derivatives  based 
on  the  geometry  specified  (Cartesian,  Cylindrical,  etc.).  This  software  tool  has  drastically 
shortened  the  programming  and  debugging  phases  of  code  development,  and  after 
reaching  its  current  mature  state  it  has  allowed  the  realization  of  a  relatively  fast 
turnaround  on  the  study  of  new  problems/problem  formulations.  In  principle,  this  code 
generation  script  could  be  turned  into  a  Mathematica  package  and  made  available  to  any 
researcher  who  employs  finite  difference  methods  to  solve  partial  differential  equations. 
As  noted  in  the  Overview  section  of  this  report,  a  partial  listing  of  the  script  is  attached  as 
an  Appendix  to  this  report. 
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Appendix  I: 

Mathematica-based  Code  Generation  Tool 

(partial  code  listing) 


Formulation-P_  Time-dependentPrimitive Variables,  nb 


Primitive  variable  formulation  of  the 
PDEs 

governing  time-dependent  binary  mixing 


Richard  Dobbins  -  January  2009 


Session 

<>  Working  directory 

Directory [ ] 

SetDirec tory [ MC : \\Documents  and  Se ttings\\rd\\My  Documents 

\\Work  -  EngineeringW-  RESEARCH  -\\PROJECTS\\Fluids\\MixingPipe" ] 

SetDirec tory [ "F : \\PRO JECTS\\Combustion\\OneStep_Dif fusion" ] 

F : \PROJECTS\Combustion\OneStep_Dif fusion 

<>  Miscellaneous 

<<  Format  .m 

Make  sure  that  the  package  FortranDf ormat . m  is  not  already  loaded,  or  else 

Fort ranAssign' s  formatting  of  some  numbers  as  DP  constants  will  not  work  correctly! 

Off [General : : " spell 1" ] 

Of f [General: : "spell” ] 

NormalPageWidth =  PageWidth  /  .  Options [$Output,  PageWidth]; 

NotZeroQ [X _ ]  ;  =  Not [Developer "ZeroQ[X]  ]  ; 

O  Error  messages 

problemsetup: ;notcomplete  = 

"Basic  arrays  (e.g.  variables,  derivatives, 

assumptionsRules)  are  not  yet  defined.  Must  complete  problem  setup."; 


<>  Page  Width  and  linebreaking 
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Definitions  &  initializations 


■  Problem  setup 


<>  Coordinate  system 

<<  Calculus'VectorAnalysis' 

Se tCoordinates [Cylindrical [r,  6,  z] ] 


Using  NEW  VectorAnalysis  package  which  handles  tensor  operations 


Cylindrical [r#  0,  z] 

<>  Variables  &  material  properties 

pres  =  P  [ t,  r ,  6,  z]  ; 
vr  =  U [ t ,  r ,  Q ,  z  ]  ; 
v0  =  V  [  t ,  r ,  Q,  z  ]  ; 
vz  =  W[t,  r ,  e ,  z ]  ; 
v  =  (vr,  v0,  vz }  ; 
temp  =  T  [  t ,  r,  6,  z]  ; 
yk  =  YK  [ t ,  r,  Q,  z]  ; 

p  =  RHO[pres,  temp,  yk]  ; 

cp  =  CP  [temp]  ; 

Acp =  LACP [ temp] ;  (*  A/cp  -  a  simple  function  of  T  ★) 

iu  =  PR  Acp ;  (*  PR  -  Prandtl  number  *) 

1 

pDk  =  - Acp;  (*  LEk  -  Lewis  number  of  species  k  *) 

LEk 

vCr  =  UC[t,  r,  0,  z]  ; 
vce  =  VC[t,  r,  6,  z]  ; 
vCz  =  WC  [t,  r,  Q,  z]  ; 
vC  =  {vCr,  vC0,  vCz }  ; 

wdot  =  SOURCEftemp,  yk]  ; 
qdot  =  HEAT  [temp,  yk]  ; 


gr  =  GR; 
g0 =  GTHETA; 
gz  =  GZ ; 

g  =  {gr ,  g0,  gz}  ; 
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O  Governing  equations 

n  =  a  (Grad[v]  +  Transpose [Grad [v] ] )  +  |/jB  -  —  Div[v]  IdentityMatrix[3  ]  ; 

(★★  deviatoric  stress  tensor  **) 

Continuity-Equation  =  dtp  +  Div[p  v]  ; 

NavierStokesEquation  =  p  at  v  +  p  v.  Grad  [v]  +  Grad[pres]  -  pg  -  Div[Il]  ; 
Energy-Equation  - 

Acp  qdot 

p  temp  +  p  v  .Grad  [temp]  -  Div  [Acp  Grad  [temp]  ] - Grad  [cp]  .Grad  [temp] - ; 

cp  cp 

SpeciesConservation  =  p  atyk  +  p  O.Grad[yk]  +  Div[pyk  vC]  -  Div [pDk  Grad  [yk]  ]  -  cjdot; 


<>  Boundary  conditions 

(*  INLET  *) 

Blp  =  NavierStokesEquation[ [3]  ]  ; 

Blu  =  vr; 

Blw  =  vz  -  Winlet; 

Bit  =  temp  -  Tinlet; 

Bly  =  yk  -  YKinlet; 

(*  SYMMETRY  *) 

B2p  =  drpres; 

B2u  =  vr; 

B2w  =  dT  vz  ; 

B2 1  =  dT temp; 

B2y  =  ar yk; 

(*  WALL  *) 

B3p  =  NavierStokesEquation[ [1]  ]  ; 

B3u  =  vr; 

B3w  =  vz; 

B3t  =  temp  -  Twall; 

B3y =  -pDk  Grad [yk]  [  [1] ]  ; 

(*  no  diffusion  of  species  k  into  the  wall:  Vk>rYk  =  0,  using  Fick'  s  Law  *) 
(*  OUTLET  ★) 

B4p  =  pres  -  AtmPressure; 

B4u=dzvr  +  vr;  (*  Robin  condition  *) 

B4w  =  dz  vz ; 

B4 1  =  dz temp; 

B4y  =  azyk; 


<>  " Rules  "  expressing  physical  assumptions 
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O  Assumptions  made  in  this  problem 

*  2-D  axisymmetnc  flow  (no  swirl) 

*  negligible  bulk  viscosity 

*  axial  gravity 

assumptionsRules  = 

{ 

(*  steadyFlowRule,  *) 
twoDimensionalAxisymmetricFlowRule, 
negligibleBulkViscosityRule, 
axialGravityRules 
}  //  Flatten; 

<>  Problem  summary 

problems  "Time-dependent  Axisymmetric  Diffusion  Flame  with  One-Step  Chemistry"; 
equations =  {"Continuity",  "Radial  Momentum", 

"Axial  Momentum",  "Temperature",  "Species  Conservation"}; 
equationnumbers  =  Range [Length[equations]  ]  ; 

boundaries =  {"Inlet",  "Axis  of  Symmetry",  "Wall",  "Outflow"}; 

variables  =  {P,  U,  W,  T,  YK}  ; 

species =  {"CH4",  "02",  "C02",  "H20",  "N2"}; 

nspecies  =  Length [species]  ;  (*  number  of  species  YK  *) 

derivtypes =  {r,  z,  rr,  zz,  rz}; 

nderivtypes  =  Length [derivtypes]  ; 

Outer [ StringJoin,  Map [ToString,  variables],  Map [ToString,  derivtypes]  ]  //  Flatten; 
derivatives  =  Map [ToExpression,  %]  ; 

Print [ 

"PROBLEM  SUMMARY;  ",  problem,  "\n", 

"Equations  =  ",  equations,  "\n", 

"Variables  =  ",  variables.  If [nspecies >  0 ,  "\n",  ""], 

If  [nspecies  >  0 ,  "Species  =  ",  ""],  If  [nspecies  >  0 ,  species,  ""] 

3; 

PROBLEM  SUMMARY;  Time-dependent  Axisymmetric  Diffusion  Flame  with  One-Step  Chemistry 
Equations  {Continuity,  Radial  Momentum,  Axial  Momentum,  Temperature,  Species  Conservation 
Variables  {P,  U,  W,  T,  YK} 

Species  { CH4  ,  02,  C02,  H20,  N2  } 
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■  Compact  Scheme  code  generation  tools 


o  Translation  Jrom  Mathematica  expressions  to  code  symbols 


derivativeTranslationRules  = 

{ 


Derivative [ 0 , 
Derivative [ 0 , 
Derivative [0, 
Derivative [ 0 , 
Derivative [ 0 , 
Derivative [ 0 , 
Derivative [ 0 , 
Derivative [0, 
Derivative [ 0 , 
Derivative [ 1, 


1,  0 ,  0]  [A_]  [t,  r ,  0,  z] 

0,  1,  0]  [A_]  [t,  r ,  0,  z] 

0,  0 ,  1]  [A_]  [t,  r,  0,  z] 

2,0,0]  [A_]  [t,  r,  0,  z] 
0,  2,  0]  [ A_]  [t,  r,  0,  z] 

0,  0,  2]  [A_]  [t,  r,  0,  z] 

1,  0,  1]  [A_]  [t,  r,  0,  z] 

1,  1,  0]  [A_J  [t,  r,  0,  z] 

0,  1,  1]  [A  ]  [t,  r,  0,  z] 

0  /  o,  0  ]  [  A_  ]  [  t ,  r ,  0 ,  z  ] 


Symbol  [ToString  [A]  <> 
Symbol  [ToString  [A]  <> 
Symbol  [ToString  [A]  <> 
Symbol  [ToString  [A]  <> 
Symbol  [ToString [A]  <> 
Symbol  [ToString  [A]  <> 
Symbol  [ToString  [A]  <> 
Symbol  [ToString  [A]  <> 
Symbol  [ToString  [A]  <> 
Symbol  [ToString  [A]  <> 


(*  "p"  =  derivative  w.r.t.  p  *) 

Derivative [  1 ,  0,  0]  [ A _ ]  [P[t,  r,  0,  z]  ,  T[t,  r,  0,  z]  ,  YK[t,  r,  6 

Symbol  [ToString  [A]  <>  npn].  Derivative  [  0 ,  1,  0]  [A_]  [P[t,  r,  6 

T  [  t ,  r,  0,  z]  ,  YK  [  t ,  r,  0,  z]  ]  Symbol  [ToString  [A]  <>  "t"]  , 

Derivative  [0,  0,  1]  [A_]  [P[t,  r,  0,  z]  ,  T[t,  r,  0,  z],  YK  [  t ,  r,  6 

Symbol [ToString [A]  <>  "y" ]  , 

Derivative  [2 ,  0,  0]  [A_]  [P[t,  r,  0,  z]  ,  T[t,  r,  0,  z]  ,  YK[t,  r,  e 
Symbol [ToString [A]  <>  "pp"] , 

Derivative [ 0 ,  2,  0]  [A  ]  [P[t,  r,  0,  z]  ,  T[t,  r,  0,  z]  ,  YK[t,  r,  e 
Symbol [ToString [A]  <>  "tt"]  , 

Derivative  [0,  0,  2]  [A_]  [P[t,  r,  0,  z]  ,  T[t,  r,  0,  z]  ,  YK  [  t ,  r,  6 
Symbol [ToString [A]  <>  "yy"] , 

Derivative [1,  1,  0]  [A_]  [P[t,  r,  0,  z]  ,  T[t,  r,  0,  z]  ,  YK [ t ,  r,  6 
Symbol [ToString [A]  <>  "pt"] , 

Derivative [0,  1,  1]  [A_]  [P[t,  r,  0,  z]  ,  T[t,  r,  0,  z]  ,  YK [ t ,  r,  6 
Symbol [ToString [A]  <>  "ty"], 

Derivative [  1,  0,  1]  [A_]  [P[t,  r,  0,  z]  ,  T[t,  r,  0,  z]  ,  YK[t,  r,  6 
Symbol [ToString [A]  <>  "py"]. 


Derivative  [  1]  [A_]  [T  [t,  r,  0,  z]]  Symbol  [ToString  [A]  <>  "t"]  , 
Derivative  [  2  ]  [A_]  [T[t,  r,  0,  z]  ]  Symbol  [ToString  [A]  <>  "tt"] 


}; 


toCodeNotation[pde_]  := 

( 

( 

Expand [pde]  /.  derivativeTranslationRules 
)  /.  { A_  [  t ,  r,  0,  z]  -♦  A} 

)  /.  { A_  [  P ,  T,  YK]  -*  A,  A_[T]  ->  A}  /.  {A_[T,  YK]  -♦A} 


"r"]  , 

«q"]  # 

»  z  "  ]  , 

"  r  r  "  ]  , 

" qq" ] , 

"  z  z  "  ]  , 
"rz"] , 
"rq"] , 
»qz" ] , 
”t"]  , 


,  z  ]  ] 
'/  z]  , 

,  Z]] 
,  Z]] 
,  Z]] 
,  Z]] 

,  z]] 
,  z]  ] 
,  z]] 
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O  Array  pointers 

(*  "pad"  pointers  with  spaces  where  needed  to  improve  code  legibility  *) 

padPointers [ptrs_]  := 

Module [ {ptrnames,  newptrnames =  {},  ptrvalues,  nptrs,  paddedlength,  i,  ptrinchars}, 
ptrnames  = ptrs [ [All ,  1]  ] / 
ptrvalues  =  ptrs [ [All,  2]  ]  ; 
nptrs  =  Length [ptrnames]  ; 

paddedlength  =  Max  [StringLength /@  ptrnames]  ; 

For[i  =  1,  i  s  nptrs,  i+  +  , 
ptrinchars  =  Characters [ptrnames [ [i] ] ] ; 

If [paddedlength  >  Length [ptrinchars] , 
ptrinchars  =  PadRight [ptrinchars,  paddedlength,  "  "  ]  ,  ] ; 

AppendTo [newptrnames,  StringJoin [ptrinchars]  ]  ; 

]  ; 

Return [Table [ {newptrnames [ [i] ] ,  ptrvalues [ [i] ]} ,  {i,  nptrs}]]; 

]; 


Module [ 

{ 

i, 

indexHeadEqns  =  "Meqn", 
indexHeadVars =  "Mvar", 
indexHeadDers =  "Kder", 
indncune, 
indindx 

If [Not [VectorQ [eguationn umbers] ]  |  | 

Not [VectorQ [variables] ]  |  |  Not [VectorQ [derivatives] ] , 

Message [problemse tup: : notcomplete]  ;  Abort [];,  ]; 

Clear [ indicesEqns ,  indicesVars,  indicesDers]  ; 

indicesEqns  =  StringJoin [ indexHeadEqns,  #]  &  /@  {ToString  /©variables) ; 
indicesVars  =  StringJoin[ indexHeadVars,  #]  &  /©  (ToString /@  variables}  ; 
indicesDers  =  StringJoin [ indexHeadDers,  #]  & /©  (ToString /@  derivatives} ; 

indicesEqns  =  Table [ (indicesEqns [ [i]  ]  ,  i},  {i.  Length [indicesEqns] }] ; 
indicesVars  =  Table [ (indicesVars [ [i]  ]  ,  i},  (i.  Length [indicesVars] }] ; 
indicesDers =  Table [ (indicesDers [ [i] ] ,  i},  {i,  Length [indicesDers] }] ; 

If [nspecies  >  1 , 

(indname,  indindx}  =  Position[indicesEqns,  "MeqnYK"]  //  Flatten; 
indicesEqns  =  ReplacePart[ 
indicesEqns, 

StringJoin [ H (/  ( " , 

ToString [indicesEqns [ [indname,  indindx + 1]  ]]  ,  "-l+k,  k=l,NSP  P}  /}" ], 


7 


Formulation-P  Time-dependentPrimitive  Variables,  nb 


{indname,  indindx+1) 

] ; 

{indname,  indindx}  =  Position[indicesVars,  "MvarYK" ]  //  Flatten; 
indicesVars  =  ReplacePart[ 
indicesVars, 

StringJoin[ " (/  (", 

ToString [indicesVars [ [indname,  indindx+1]]],  "-1+k,  k=l,NSP_P)  /)* ], 
{indname,  indindx+1) 

]; 

YKdernames  =  Select [indicesDers [ [All,  1]],  StringMatchQ [ # ,  "*YK*"]  &]  ; 
indlist  =  Flatten[Posi tion[ indicesDers,  #]  & /@  YKdernames,  1]; 

Do  [ 

{indname,  indindx)  = indlist [ [k] ] ; 
indicesDers  =  ReplacePart[ 
indicesDers, 

StringJoin [ " (/  ( " ,  ToString [ indicesDers [ [indname,  indindx+1]]], 

"+ (k-1) *NDRTP_P,  k=l,NSP_P)  /)"], 

{indname,  indindx+1) 

]  /*  / 

{k.  Length [indlist] ) 

] ; 

,  (*  ELSE  *) 

] ; 


indicesEqns  =  indicesEqns  //  padPointers; 
indicesVars  =  indicesVars  //  padPointers; 
indicesDers  =  indicesDers  //  padPointers; 


]  ; 


O  Code  symbols 

Options [wri teCodeSymbols]  = 

{ 

codeform  -+  FortranForm, 
indent  -»  "  "  , 

tocode  -+  {  ) 

) ; 


wri teCodeSymbols: susage  = 

MwriteCodeSymbols[ symbol array, output, options  ] : \n 

Writes  the  symbols  (variables,  derivatives,  etc.)  in  a  given 

PDE  in  code  form.  Options  work  as  in  Jacobian  components  routines."; 
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writeCodeSymbols [symbolarray^,  output_,  options _ ]  := 

Module [ 

{ 

lhs,  rhs,  eqn,  sublistlength, 

sublistwidth,  symbol,  symboll,  symbol2,  symbol3,  symbol inchars, 
paddedlength  =  4 ,  (★★  number  of  chars/spaces  between  1st  char  and  "  ="  ★*) 

MyCodeForm = codeform  /.  {options}  /.  Options [wri teCodeSymbols] , 
MyCodelndent = indent  /.  {options}  /.  Options [writeCodeSymbols] , 

MyToCode =  tocode  /.  {options}  /.  Options [writeCodeSymbols] 

}/ 

sublistwidth  =  Dimensions [symbolarray]  [ [1] ] ; 

If  [sublistwidth  ==  1, 
symbol  =  ToString [symbolarray [ [1]  ]  ]  ; 

If  [output  ==  1  |  |  output  ==  3  ,  Print  [symbol,  "  =  "  ,  "  \n"],]; 
symbolinchars  =  Characters [ symbol]  ; 

If [Length [symbolinchars]  <  paddedlength, 
symbol  =  StringJoin [PadRight [symbolinchars,  paddedlength,  "  "]],]; 
lhs  =  StringJoin [MyCodelndent,  symbol,  "  =  "]; 

MyToCode  =  Flatten [Append [MyToCode,  {"\n",  lhs}] ] ; 

,  (★*  ELSE  ★★) 

symboll =  ToString [symbolarray [ [1] ] ] ; 
symbol2  =  ToString [symbolarray [ [2]  ]  ]  ; 
symbol3  =  ToString  [symbolarray  [  [3]  ]  ]  ,* 

If  [output  ==  1  |  |  output  ==  3  , 

Print  [symboll,  "  =  "  ,  symbol2,  "  (",  symbol3,  "  ,  ind)  "  ,  "\nM  ],]; 
symbolinchars  =  Characters [ symboll]  ; 

If [Length [symbolinchars]  <  paddedlength, 
symboll  =  StringJoin [PadRight [symbolinchars,  paddedlength,  "  "]],]; 
lhs  =  StringJoin[MyCodeIndent,  symboll,  "  =  "]; 
rhs = StringJoin [symbol2,  "(",  symbol3,  ",ind)"]; 

MyToCode =  Flatten [Append [MyToCode,  {"\nn,  lhs,  rhs}]]; 

] ; 


If [FreeQ [ {options } ,  tocode]. 

If  [  (output  ==  2  |  |  output  ==  3)  ,  Print  @@  MyToCode,  ]  , 
Re  turn [ MyToCode ] 

]  ; 


]; 
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<>  Equation  residuals 

Options [writeEquationResidual]  = 

{ 

codeform->  FortranAssign, 
indent  -»  "  " , 

tocode  -»  { } 

}; 


writeEquationResidual: : usage  = 

"writeEquationResidual [pde , pdeindex, output, options _ ] : \n 

Writes  the  equation  residual  for  a  given  PDE 

in  code  form.  Options  work  as  in  Jacobian  components  routines."; 

wri teEquationResidual [pdein_,  pdeindex_,  output,  options _ ]  := 

Module [ 

{ 

lhs,  rhs,  eqn, 

MyCodeForm =  codeform  /.  {options}  /.  Options [writeEquationResidual]  , 
MyCodelndent =  indent  /.  {options}  /.  Options [writeCodeSymbols] , 

MyToCode =  tocode  /.  {options}  /.  Options [wri teEquationResidual] 

}, 

pde  =  pdein; 

If  [output  s*  1  |  |  output  ss  3,  Print  ["F  =  " ,  pde,  "\n"],]; 

lhs  = StringJoinfMyCodelndent,  "EQ0(",  indicesEqns [ [pdeindex,  1]],  ",ind)  =  "]; 
rhs  =  MyCodeForm [pde, 

AssignBreak ->  {1000,  "\n  &"}/ 

Assignlndent  ->  99 , 

AssignOptimize  ->  False, 

AssignPrecision ->  Infinity, 

AssignTemporary ->  {"tmpO" ,  Sequence} 

][[i]]  ; 

MyToCode =  Flatten [Append [MyToCode,  {"\n",  lhs,  rhs}]]; 

If [FreeQ [ {options } ,  tocode], 

If  [  (output  ==  2  |  |  output  ==  3)  ,  Print  @@  MyToCode,  ]  , 

Return [MyToCode] 

]; 


]; 
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O  Local  components  of  the  Jacobian 

Options  [  j  acobianLocalComponents]  = 

{ 

codeform-*  FortranAssign, 
indent  -»  ■  ", 

tocode  ->  { } 

}; 


j acobianLocalComponents: :usage  = 

"j  acobianLocalComponents [pde, pdeindex, vars , output, options _ ] : \n 

★  'pde'  is  the  equation  in  symbolic  form;  1 

pdeindex'  is  the  number  of  this  PDE  in  the  system  of  PDEs\n 

★  'vars'  is  a  list  of  all  possible  variables  (unknowns) 

which  can  occur  in  the  equations\n 

★  'output'  is  a  flag  which  determines  whether  Mathematica  summaries  of 

the  components  or  actual  code  is  written  (1  =  Mathematica  only, 

2  =  code  only,  3  =  both) ;  current  settings  generate  Fortran  code\n 

★  'options'  is  a  sequence  of  rules  (not  required)"; 

j acobianLocalComponents [pdein_,  pdeindex_,  vars_,  output_,  options  ]  := 

Module [ 

{ 

lhs,  rhs, 

MyCodeForm  =  codeform  /.  {options}  /.  Options [ j acobianLocalComponents]  , 
MyCodelndent =  indent  /.  {options}  /.  Options [writeCodeSymbols] , 

MyToCode =  tocode  /.  {options}  /.  Options [ j acobianLocalComponents] 

pde  =  pdein//.  j acobianldealGasRule; 

Do  [ 

var  =  vars [ [k] ] ; 
tmpC)  -  D[pde,  var]  ; 

If  [var  —  —  —  T,  tmpl  =  LACP t  D [pde ,  LACP]  ,  tmpl  =  0]  ; 

If  [var  =  =  =  T,  tmp2  =  CPt  D[pde,  CP]  ,  tmp2  =  0]  ; 

If  [var  =  =  =  T,  tmp3  =  LACPtt  D[pde,  LACP t]  ,  tmp3  =  0]  ; 

If  [var  =  -  -  T,  tmp4  =  CPtt  D  [pde ,  CPt]  ,  tmp4  =  0]  ; 

If  [var  =  =  =  YK,  tmp5  =  WMy  D[pde,  WM]  ,  tmp5  =  0]  ; 

DJ [pdeindex,  k]  = 

(  (tmpO  +  tmpl  +  tmp2  +  tmp3  +  tmp4  +  tmp5)  //  Expand)  //  .  j  acobianldealGasRes toreRules; 

If  [ 

DJ [pdeindex,  k]  =!=  0, 

If  [output  ==  1  |  |  output  ==  3 , 

Print ["d  F",  ToString [pdeindex] , 

"  /  d  (",  vars[  [k]  ]  ,  ")  =  ",  DJ  [pdeindex,  k]  ,  "\n"],]; 
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lhs =  StringJoin [MyCodelndent,  "DJ(",  indicesEqns [ [pdeindex,  1]], 
"  ,  "  ,  indicesVars [ [k,  1]],  " , ind)  =  "]; 
rhs  =  MyCodeForm[DJ [pdeindex,  k]  , 

AssignBreak  -►  {  1000  ,  M\n  &”}, 

Assignlndent -> 

AssignOptimize -►  False, 

AssignPrecision  ->  Infinity, 

AssignTemporary ->  {" tmpO",  Sequence} 

][[i]]  ; 

MyToCode =  Flatten [Append [MyToCode,  {"\n”,  lhs,  rhs}]]; 

,  (*  ELSE  *) 

]  ; 


,  {k.  Length [vars]  } 


]  ; 


If [FreeQ [ {options} ,  tocode] , 

If  [  (output  “  2  |  |  output  =*  3)  ,  Print  @@  MyToCode,  ]  , 
Re  turn [ MyToCode] 

] ; 


]; 


<>  Time  derivative  (diagonal)  terms  of  the  Jacobian 

Options  [  j  acobianTimeDerivatives]  = 

{ 

codeform  -►  FortranAssign, 
indent  *+  "  ”, 

tocode  -►  {  } 

}; 


j acobianTimeDerivatives:  : usage  = 

"j acobianTimeDerivatives [pde,pdeindex, vars, output, options _ ] : \n 

★  'pde'  is  the  equation  in  symbolic  form;  1 

pdeindex'  is  the  number  of  this  PDE  in  the  system  of  PDEs\n 

★  'vars'  is  a  list  of  all  possible  variables  (unknowns) 

which  can  occur  in  the  equations\n 

*  'output1  is  a  flag  which  determines  whether  Mathematics  summaries  of 

the  components  or  actual  code  is  written  (1  =  Mathematics  only, 

2  =  code  only,  3  =  both) ;  current  settings  generate  Fortran  code\n 

*  'options'  is  a  sequence  of  rules  (not  required)”; 

j acobianTimeDer ivatives [pdein  ,  pdeindex  ,  vars  ,  output  ,  options  _ ]  :  = 

Module [ 


lhs,  rhsl,  rhs2, 

MyCodeForm = codeform  /.  {options}  /.  Options [j acobianLocalComponents]  , 


{ 
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MyCodelndent =  indent  /.  {options}  /.  Options [writeCodeSymbols] , 
MyToCode =  tocode  /.  {options}  /.  Options [j acobianLocalComponents] , 
tders,  tmptders 


tmptders  =  Outer [  StringJoin,  Map [ToString,  vars] ,  {"t"}  ]  //  Flatten; 
tders  =  Map [ToExpression,  tmptders] ; 

Remove [ tmptders]  ; 

pde  =  pdein//.  j acobianldealGasRule; 

Do  [ 

tmpsDtpde,  tders  [[k]]]; 

DJT [pdeindex,  k]  =  tmp//.  j acobianldealGasRes toreRules; 

If  [ 

DJT [pdeindex,  k]  =1=0, 

If  [output  ==  1  |  |  output  ==  3, 

Print ["d  FM,  ToString [pdeindex] , 

"  /  d  (",  tders  [  [k]  ]  ,  ")  =  ",  DJT  [pdeindex,  k]  ,  " \n"],]; 

lhs  = StringJoin [MyCodelndent,  "DJ(",  indicesEqns [ [pdeindex,  1]], 
indice sVars [ [k,  1]],  " , ind)  =  "]; 
rhsl = StringJoin[ "DJ ( " ,  indicesEqns [ [pdeindex,  1]], 
indicesVars [ [k,  1]],  " , ind)  +  ■]; 
rhs2  s  MyCodeForm[ DJT [pdeindex,  k]  ddt, 

AssignBreak  -*  {1000,  "\n  &"}# 

Assignlndent  -* 

AssignOptimi ze  ->  False, 

AssignPrecision -►  Infinity, 

AssignTemporary {"tmpO",  Sequence} 

][[!]]  ; 

MyToCode =  Flatten [Append [MyToCode,  {"\n",  lhs,  rhsl,  rhs2}]]; 

,  (*  ELSE  ★) 

]  ; 


,  {k,  Length [vars] } 


]; 


If [FreeQ [ {options} ,  tocode]. 

If  [  (output  ==  2  |  |  output  ==  3)  ,  Print  @@  MyToCode,  ]  , 
Return [MyToCode] 

]  ; 


]; 


<>  Spatial  components  of  the  Jacobian 
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Options [ jacobianSpatialComponents]  = 

{ 

codefonn-»  FortranAssign, 
indent  ->  "  ", 

tocode  ->  {  } 

}; 


jacobianSpatialComponents: : usage  = 

" j  acobianSpatialComponents[pde, pdeindex, derivs,  output , options _ ] : \n 

★  'pde'  is  the  equation  in  symbolic  form;  1 

pdeindex'  is  the  number  of  this  PDE  in  the  system  of  PDEs\n 

★  'derivs'  is  a  list  of  all  possible  derivatives  which  can  occur  in  the  equations\n 

★  'output'  is  a  flag  which  determines  whether  Mathematica  summaries 

of  the  components  or  actual  code  is  written  (1  =  Mathematica  only, 

2  =  code  only,  3  =  both) ;  current  settings  generate  Fortran  code\n 

★  'options'  is  a  sequence  of  rules  (not  required)"; 
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j acobianSpatialComponents[pdein_,  pdeindex,  derivs_,  output_,  options _ ]  := 

Module [ 

{ 

lhs,  rhs ,  eqn, 

MyCodeForm =  codeform  /.  {options}  /.  Options [j acobianSpatialComponents] , 
MyCodelndent =  indent  /.  {options}  /.  Options [writeCode Symbol s] , 

MyToCode =  tocode  /.  {options}  /.  Options [ j acobianSpatialComponents] 

h 

pde  =  pdein//.  j acobianldealGasRule; 

Do  [ 

tmp  =  D [pde,  derivs [ [k]  ]  ]  ; 

DQ[pdeindex,  k]  =  tmp//.  j acobianldealGasRestoreRules; 

If  [ 

DQ[pdeindex,  k]  =!=  0, 

I  f  [output  ==  1  |  |  output  ==  3 , 

Print["d  F",  ToS tring[pdeindex]  , 

11  /  d  (",  derivs  [  [k]  ]  ,  ")  =  ",  DQ [pdeindex,  k]  ,  "  \n"  ],]; 

lhs =  StringJoin [MyCodelndent,  "DQ(",  indicesEqns [ [pdeindex,  1]], 

11 ,  " ,  indicesDers  [  [k,  1]],  "  ,  ind)  =  "]; 
rhs  =  MyCodeForm [DQ [pdeindex,  k]  , 

AssignBreak -►  { 1000,  "\n  &"}/ 

Assignlndent ->  11 "  , 

AssignOptimize  -►  False, 

AssignPrecision ->  Infinity, 

AssignTemporary ->  {" tmpO",  Sequence} 

][[!]]  ; 

MyToCode =  Flatten [Append [MyToCode,  {"\n",  lhs,  rhs}]]; 

,  (*  ELSE  *) 

]  ; 


,  {k,  Length[derivs] } 


]? 


If [FreeQ [ {options} ,  tocode]. 

If  [  (output  ==  2  |  |  output  ==  3 )  ,  Print  @@  MyToCode,  ]  , 
Return [MyToCode] 

] ; 


]; 
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RESIDUALS  &  JACOBIAN  -  INTERIOR  POINTS 


■  Pressure 


O  Residuals 


pde  =  Continuity-Equation; 

pde  =  pde  /  .  assumptionsRules //  Expand; 

pde  =  (r  *  pde)  //  Expand; 

PDE[1]  =  toCodeNotation[pde]  //.  equationsIdealGasRules 
Remove [pde] ; 


Pt  r  RHO 
P 

Pz  r  RHO  W 
P 


r  RHO  Tt 


T 

r  RHO  Tz  W 


RHO  U  • 


Pr  r  RHO  U  r  RHO  Tr  U 


+  r  RHO  Wz 


r  RHO  U  WM  YKr 
WK 


+  r  RHO  Ur  + 

r  RHO  WM  YKt 
WK 


r  RHO  W  WM  YKz 
WK 


O  Jacobian 


j acobianLocalComponents [PDE [1]  ,  1 ,  variables,  1]; 
j acobianTimeDerivatives [PDE [1]  ,  1,  variables,  1]; 
j acobianSpatialComponents[PDE [1]  ,  1,  derivatives,  1]  ; 


d  Fi  /  d  (P)  = 
r  RHO  Tz  W 


r  RHO  Tt  RHOU  r  RHO  Tr  U  r  RHO  Ur 


P  T  P 

r  RHO  Wz  r  RHO  U  WM  YKr 


P  T  P 

r  RHO  WM  YKt  r  RHO  W  WM  YKz 


P  T 


P  WK 


P  WK 


P  WK 


d  FI  /  d  (U)  =  RHO  - 


Pr  r  RHO  r  RHO  Tr  r  RHO  WM  YKr 


d  Fl  /  d  W)  = 


Pz  r  RHO  r  RHO  Tz  r  RHO  WM  YKz 


WK 


d  Fl  /  d  (T) 

Pz  r  RHO  W 
PT 


Pt  r  RHO  2  r  RHO  Tt  RHOU  Pr  r  RHO  U  2  r  RHO  Tr  U  r  RHO  Ur 


P  T 

2  r  RHO  Tz  W 
T2 


T2 

r  RHO  Wz 


T  P  T  T2  T 

r  RHO  U  WM  YKr  r  RHO  WM  YKt  r  RHO  W  WM  YKz 


T  WK 


T  WK 


T  WK 


d  Fl  /  d  (YK) 

Pz  r  RHO  W  WM 
P  WK 


Pt  r  RHO  WM  r  RHO  Tt  WM  RHO  U  WM  Pr  r  RHO  U  WM  r  RHO  Tr  U  WM  r  RHO  Ur  WM 
P  WK  +  T  WK  WK  P  WK  +  T  WK  WK 

r  RHO  Tz  W  WM  r  RHO  WM  Wz  2  r  RHO  U  WM2  YKr  2  r  RHO  WM2  YKt  2  r  RHO  W  WM2  YKz 
T  WK  WK  WK2  WK2  WK2 


d  Fl  /  d  (Pt) 


r  RHO 
P 


r  RHO 


d  Fl  /  d  (Tt) 


T 
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d  FI  /  d  (YKt) 

d  FI  /  d  (Pr)  - 

d  FI  /  d  (Pz)  - 

d  FI  /  d  (Ur)  - 

d  FI  /  d  (Wz)  - 

d  FI  /  d  (Tr)  - 

d  FI  /  d  (Tz)  - 

d  FI  /  d  (YKr )  = 

d  FI  /  d  ( YKz )  = 


r  RHO  WM 
WK 


r  RHOU 
P 


r  RHO  W 
P 


r  RHO 


r  RHO 


r  RHO  U 
T 


r  RHO  W 
T 


r  RHO  U  WM 
WK 


r  RHO  W  WM 
WK 


■  Radial  velocity 


<>  Residuals 

pde  =  NavierStokesEquation[  [1]  ]  ; 

pde  =  pde  /.  assumptionsRules  //  FullSimplify  //  Expand; 
pde  =  (r2  *  pde)  //  Expand; 

PDE [2]  =  toCodeNotationfpde]  //.  equationsIdealGasRules 
Remove [pde] ; 

Pr  r2  +  4  LA-P  PR  u  +  lac pt  pr  r  Tr  U  -  —  LACP  PR  r  Ur  -  —  LACPt  PR  r2  Tr  Ur  + 

3  3  3  3 

r2  RHO  U  Ur  -  y  LACP  PR  r2  Urr  +  r2  RHO  Ut  -  LACPt  PR  r2  Tz  Uz  -  LACP  PR  r2  Uzz  + 
r2  RHO  Uz  W  -  LACPt  PR  r2  Tz  Wr  -  y  LACP  PR  r2  Wrz  +  y  LACPt  PR  r2  Tr  Wz 


Formulation-P _ Time-dependen tPrim iti ve  Variables  n  b 
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<?  Jacobian 


j a c obi an Local Components [PDE [2 ] ,  2, 
jacobianTimeDerivatives[PDE [2]  ,  2, 
j acobianSpatialComponents[PDE [2]  , 


variables,  1]  ; 
variables,  1]  ; 
2,  derivatives. 


i]; 


d  F2  /  d  (P) 

d  F2  /  d  (U) 

d  F2  /  d  (W) 


r2  RHO  U  Ur 
P 


4  LACP  PR 

3 


r2  RHO  Uz 


2_ 

3 


r2  RHO  Ut  r2  RHO  Uz  W 
P  +  P 

LACPt  PR  r  Tr  +  r2  RHO  Ur 


d  F2  /  d  {T) 


4  LACPt  PR  U 


2  4 

—  LACPt  t  PR  r  Tr  U  -  y  LACPt  PR  r  Ur  - 


y  LACPt  t  PR  r2  Tr  Ur  -  r  RHO  U  Ur  _  4_  LACpt  pR  r2  Urr  _  £  RHO  Ut  _  LACptt  pR  r2 

r2  RHO  Uz  W  12 

LACPt  PR  r2  Uzz  -  - - - LACPtt  PR  r2  Tz  Wr  -  y  LACPt  PR  r2  Wrz  +  y  LACPtt 


Tz  Uz  - 
PR  r2  Tr  Wz 


d  F2  /  d  (YK)  = 


r2  RHO  UUrWM  r2  RHO  Ut  WM  r2  RHO  UzWWM 


WK 


WK 


WK 


d  F2  /  d  (Ut)  =  r2  RHO 


d  F2  /  d  (Pr)  =  r2 


d  F2  /  d  (Ur) 


4  4 

y  LACP  PR  r  -  y  LACPt  PR  r2  Tr  +  r2  RHO  U 


d  F2  Id  (Uz)  -LACPt  PR  r2  Tz  +  r2  RHO  W 


d  F2  /  d  (Urr)  = 

d  F2  /  d  (Uzz)  = 

d  F2  /  d  (Wr)  - 

d  F2  /  d  (Wz)  = 

d  F2  /  d  (Wrz) 

d  F2  /  d  (Tr)  - 


-  y  LACP  PR  r2 

-LACP  PR  r2 

-LACPt  PR  r2  Tz 

2  , 
y  LACPt  PR  r2  Tr 

-  y  LACP  PR  r2 

y  LACPt  PR  r  U  -  y  LACPt  PR  r2  Ur  +  y  LACPt  PR  r2  Wz 


d  F2  /  d  (Tz) 


-  LACPt  PR  r2  Uz  -  LACPt  PR  r2  Wr 


Formulation-P  T ime-dependentPrim  itive  V anables.  nb 


18 


■  Axial  velocity 

<>  Residuals 

pde  =  NavierStokesEquation[ [3]  ]  ; 

pde  =  pde  /  .  assumptionsRules  //  FullSimplify  //  Expand; 
pde  =  (r  *  pde)  //  Expand; 

PDE [ 3 ]  =  toCodeNotationfpde]  //.  equationsIdealGasRules 
Remove [pde] ; 

Pz  r  -  GZ  r  RHO  +  ~  LACPt  PR  Tz  U  +  y  LACPt  PR  r  Tz  Ur  -  y  LACP  PR  r  Urz  - 


LACP  PR  Uz 
3 


-  LACPt  PR  r  Tr  Uz  -  LACP  PR  Wr  -  LACPt  PRrTr  Wr  +  r  RHO  UWr- 


4 


4 


LACP  PR  r  Wrr  +  r  RHO  Wt - LACPt  PRrTzWz  +  r  RHO  WWz - LACP  PRrWzz 


3 


3 
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O  Jacobian 

jacobianLocalComponents[PDE [3 ]  ,  3,  variables,  1]  ; 
j acobianTimeDerivatives [PDE [3 ]  ,  3,  variables,  1]; 
jacobianSpatialComponents[PDE[3]  ,  3,  derivatives,  1] 


d  F3  /  d  (P)  = 

GZ  r  RHO  r  RHO  U  Wr  r  RHO  Wt  r  RHO  W  Wz 

- +  ■  +  -  +  ■ 

P  P  P  P 

d  F3  /  d  (U)  = 

2  LACPt  PR  Tz  ntt^TT 

- - -  +  r  RHO  Wr 

d  F3  /  d  (W)  - 

r  RHO  Wz 

d  F3  /  d  (T)  - 

GZ  r  RHO  2  2 

- - -  +  y  LACPtt  PR  Tz  U  +  y  LACPtt  PR  r  Tz  Ur  - 

y  LACPt  PR  r  Urz  -  —————————  -  LACPtt  PR  rTrUz-  LACPt  PR  Wr  -  LACPtt  PR  rTrWr- 


d  F3  /  d  (YK)  - 

GZ  r  RHO  WM  r  RHO  UWMWr  r  RHO  WM  Wt  r  RHO  W  WM  Wz 

WK  WK  WK  WK 

d  F3  /  d  (Wt)  = 

r  RHO 

d  F3  /  d  (Pz)  -- 

r 

d  F3  /  d  (Ur) 

y  LACPt  PR  r  Tz 

d  F3  /  d  (Uz) 

LACP  PR  T  _  _  m 
- - - LACPt  PR  r  Tr 

d  F3  /  d  (Urz) 

:  -  y  LACP  PR  r 

d  F3  /  d  (Wr)  = 

-  LACP  PR  -  LACPt  PR  r  Tr  +  r  RHO  U 

d  F3  /  d  (Wz)  = 

4 

-  y  LACPt  PR  r  Tz  ♦  r  RHO  W 

d  F3  /  d  (Wrr ) 

-  LACP  PR  r 

d  F3  /  d  (Wzz ) 

=  -  y  LACP  PR  r 

d  F3  /  d  (Tr)  = 

-  LACPt  PR  r  Uz  -  LACPt  PR  r  Wr 

d  F3  /  d  (Tz)  = 

2  LACPt  PR  U  2  4 

- - -  +  y  LACPt  PR  r  Ur  -  y  LACPt  PR  r  Wz 

3 
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■  Temperature 

c>  Residuals 

pde  =  EnergyEquation; 

pde  =  pde  /  .  assumptionsRules  //  Expand; 
pde  =  (r  *  pde)  //  Expand; 

PDE [4]  =  toCodeNotation[pde]  //.  equationsIdealGasRules 
Remove [pde] ; 


HEAT  r 
CP 


-  LACP  Tr  - 


CPt  LACP  r  Tr2 
CP 


-  LACPt  r  Tr2  -  LACP  r  Trr  + 


r  RHO  Tt  - 


CPt  LACP  r  Tz2 
CP 


-  LACPt  r  Tz2  -  LACP  r  Tzz  +  r  RHO  TrU  +  r  RHO  Tz  W 


21 
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O  Jacobian 


j acobianLocalComponents [PDE [4 ]  ,  4,  variables,  1]; 
j acobianTimeDerivatives [PDE [4 ]  ,  4,  variables,  1]; 
jacobianSpatialComponents[PDE [4]  ,  4,  derivatives,  1]  ; 

r  RHO  Tt  r  RHO  TrU  r  RHO  Tz  W 

-  +  -  +  - 

P  P  P 

r  RHOTr 


d  F4  /  d  (P)  - 

d  F4  /  d  (U) 


d  F4  /  d  (W)  =  r  RHO  Tz 


d  F4  /  d  (T)  - 

CPt  LACPt  r  Tr2 
CP 


CPt  HEAT  r  m  CPt2  LACP  r  Tr2  CPtt  LACP  r  Tr2 

-  LACPt  Tr  +  - 


CP2 


LACPt t  r  Tr  -  LACPt  r  Trr 


CP2  CP 

r  RHO  Tt  CPt2  LACP  r  Tz2 


CPtt  LACP  r  Tz2  CPt  LACPt  r  Tz2 


CP 


d  F4  /  d  (YK) 


CP 


LACPt t r  Tz  -  LACPt  r  Tzz  • 


CP2 

r  RHOTr  U 


r  RHO  Tt  WM  r  RHO  Tr  U  WM  r  RHO  Tz  W  WM 


WK 


d  F4  /  d  (Tt)  -  r  RHO 


2  CPt  LACP  r  Tr 

d  F4  /  d  (Tr)  -  LACP  -  - - -  2  LACPt  r  Tr  +  r  RHO  U 


2  CPt  T.ACP  r  T7 

d  F4  /  d  (Tz)  -  cp  ~  2  ^CPt  r  Tz  +  r  RHOW 


d  F4  /  d  (Trr)  - LACP  r 


r  RHO  Tz  W 
T 


d  F4  /  d  (Tzz) 


-LACP  r 
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■  Species 


c>  Residuals 


pde  =  SpeciesConservation; 

pde  =  pde  /  .  assumptionsRules  //  Expand; 

pde  =  (r  *  pde)  //  Expand; 

PDE [5]  =  toCodeNotation[pde]  //.  equationaldealGasRules 
Remove [pde] ; 


-  r  SOURCE  +  RHO  UC  YK  + 


Pr  r  RHQ  UC  YK  r  RHO  Tr  UC  YK 


r  RHO  UCr  YK  - 


Pz  r  RHO  WC  YK  r  RHO  Tz  WC  YK 


P  T 

r  RHO  U  YKr  +  r  RHO  UC  YKr 
LACPt  r  Tz  YKz 


r  RHO  WCz  YK- 


LACP  YKr  LACPt  r  Tr  YKr 


r  RHO  UC  WM  YK  YKr 


LEk 

LACP  r  YKrr 


WK 


LEk 


+  r  RHO  W  YKz  +  r  RHO  WC  YKz  + 


LEk 

r  RHO  WC  WM  YK  YKz 


LEk 

r  RHO  YKt  - 

LACP  r  YKzz 


WK 


LEk 


<>  Jacobian 


jacobianLocalComponents [PDE [5] /  5,  variables,  1]; 
jacobianTimeDerivatives [PDE [5]  ,  5,  variables,  1]  ; 
j acobianSpatialComponents [PDE [5]  ,  5,  derivatives,  1]  ; 


d  F5  /  d  (P)  = 

r  RHO  UC  YKr 
P 


RHO  UC  YK  r  RHO  Tr  UC  YK  r  RHO  UCr  YK 
P  PT  +  P 

r  RHO  UC  WM  YK  YKr  r  RHO  YKt  r  RHO  W  YKz 
P  WK  f  P  +  P 


r  RHO  Tz  WC  YK 
PT 

r  RHO  WC  YKz 
+  — . — _ — - 

P 


r  RHO  WCz  YK  r  RHO  U  YKr 

p  +  'p" 

r  RHO  WC  WM  YK  YKz 
P  WK 


d  F5  /  d  (U)  =  r  RHO  YKr 


d  F5  /  d  (W)  =  r  RHO  YKz 


F5  /  d  (T)  = 

Pz  r  RHO  WC  YK 
P~T 

r  RHO  U  YKr 


RHO  UC  YK  Pr  r  RHO  UC  YK  2  r  RHO  Tr  UC  YK  r  RHO  UCr  YK 


2  r  RHO  Tz  WC  YK 


P  T 

r  RHO  WCz  YK 


T2 

r  RHO  UC  YKr 


r  RHOUC  WM  YK  YKr 


LACPt  t  r  Tz  YKz  r  RHO  W  YKz 


T  WK 

r  RHO  WC  YKz 


T2 

LACPt  YKr 
LEk 

LACPt  r  YKrr 
LEk 

r  RHO  WC  WM  YK  YKz 


LEk 


LACPt  t  r  Tr  YKr 
LEk 

r  RHO  YKt 


LACPt  r  YKzz 
LEk 


d  FB  /  d  (YK)  = 
r  RHO  TZ  WC 


RHO  UC 


Pr  r  RHO  UC  r  RHO  Tr  UC 


r  RHO  UCr 


Pz  r  RHO  WC 


r  RHO  WCz 


RHO  UC  WM  YK  Pr  r  RHO  UC  WM  YK  r  RHO  Tr  UC  WM  YK 


r  RHO  UCr  WM  YK 


WK 

PZ  r  RHO  WC  WM  YK 


WK 

2  r  RHO  UC  WM2  YK  YKr 


P  WK 

r  RHO  WM  YKt 


P  WK 

r  RHO  Tz  WC  WM  YK 
T  WK 

r  RHO  W  WM  YKz 


T  WK 

r  RHO  WCz  WM  YK 
WK 


r  RHO  U  WM  YKr 


WK 


2  r  RHO  WC  WM2  YK  YKz 


WK 


d  F5  /  d  (YKt) 


r  RHO 
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d  F5  /  d  (Pr) 


d  F5  /  d  (Pz) 


d  F5  /  d  (Tr) 


d  F5  /  d  (Tz) 


r  RHO  UC  YK 
P 


r  RHO  WC  YK 
P 


i  RHO  UC  YK  LACPt  r  YKr 
T  LEk 


r  RHO  WC  YK  LACPt  r  YKz 
T  LEk 


d  F5  /  d  (YKr) 


LACP 

LEk 


LACPt  r  Tr 
LEk 


+  r  RHO  U  +  r  RHO  UC  + 


r  RHO  UC  WM  YK 
WK 


d  F5  /  d  (YKz) 


LACPt  r  Tz 
LEk 


+  r  RHO  W  +  r  RHO  WC  + 


r  RHO  WC  WM  YK 
WK 


d  F5  /  d  (YKrr ) 


LACP  r 
LEk 


d  F5  /  d  (YKzz) 


RESIDUALS  &  JACOBIAN  -  INLET  BOUNDARY 


RESIDUALS  &  JACOBIAN  -  AXIS  OF  SYMMETRY 


RESIDUALS  &  JACOBIAN  -  WALL  BOUNDARY 


RESIDUALS  &  JACOBIAN  -  OUTLET  BOUNDARY 


Output  code  material  for  point.f90 


Output  code  material  for  orj.f 


Appendix  II: 

Fortran  Code  for 

RESIDUALS  and  JACOBIAN  Subroutines 

of  the  Implicit-Compact  Solver 

(partial  listing  of  a  code  created  using  the 
software  tool  presented  in  Appendix  1) 


! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


1  %%  %% 

1%%  %% 

!%%  RESIDUALS  %% 

!  %%  %% 


! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SUBROUTINE  RESIDUALS  (ISTEP) 


USE  DimensioningParameters 
USE  ProblemParameters 

USE  MethodParameters ,  ONLY  :  iDSPACE, iDTIME,  & 

Sc  TFLOW ,  DTO  ,  DTI ,  DT2  ,  Sc 

Sc  KVFIX,  NFIX 

USE  Pointers 

USE  DiscretizedProblem  ,  ONLY  :  NVAR ,  NSP,  NR ,  NZ,  NODES ,  NEL,  Sc 

Sc  RI  |  Z  J  t  Sc 

Sc  SF,  S,  SI,  S2,  SX,  Sc 

Sc  SD ,  Sc 

Sc  EQO ,  EQ1 ,  CT ,  F,  SB1,  Sc 

Sc  VELinlet 

USE  Differentiation 
USE  OneStepChemistry_Methane 
USE  Simplif iedTransport 


USE  Time_IO_Debug,  ONLY  :  WTRES 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 


double  precision  dSdt (NVAR_P/ NODES_P) 

double  precision  YK(NSP_P) ,  YKt(NSPJP),  YKinlet (NSP_P) 

double  precision  YKr(NSP_P)/  YKz  (NSP_P)  ,  YKrr(NSP_P)/  YKzz(NSPJ) 

double  precision  LACP,  LACPt ,  LEK(NSP_P)  #  DK(NSP_P)/  DKt(NSP_P)/  Sc 

!  stoichiometric  coeffs  (#  mols  produced/destroyed) 

Sc  NUK  (NSP_P)  ,  Sc 

\  mass  coeffs:  NUK(k)*WK(k) 

Sc  NUKWK(NSP_P)  ,  Sc 

Sc  WK  ( NS  P_P )  ,  SOURCE  ( NS  P_P ) 

save  isavsd  !  should  be  initialized  to  zero  by  compiler 
save  igetrsf  !  should  be  initialized  to  zero  by  compiler 
save  ioutflag 

! . upwinding . 

double  precision  CONVRDER (NVAR_P , NODE S_P) , CONVZDER (NVAR_P , NODE S_P) 
double  precision  CONVECT (NVAR_P) 

! . upwinding . 

iUPWIND  =1  !  1  =  upwind;  2  =  2nd  order  centered  (convective  terms) 


timephase  =  2*PI *TFLOW/Period 


INITIALIZATION 


call  DZERO ( F , NEL) 
call  DZERO (EQO , NEL) 
call  DZERO ( dSdt , NEL ) 
call  DZERO (CT, NEL) 


upwinding . 

call  DZERO (CONVRDER, NEL) 
call  DZERO (CONVZDER, NEL) 
call  DZERO ( CONVECT , NVAR ) 

if  (iSTEP.eq.O)  then 
call  DZERO (EQ1, NEL) 
endif 


SPATIAL  DERIVATIVES 


if  ( iDSPACE . eq . 1 )  then 
call  SPATI AL_LO  (S) 
else 

call  SPATIAL_CS  (S) 
endif 


RESIDUAL  -  INTERIOR  NODES 


DO  j  =  2,  NZ-1 
DO  i  =  2 ,  NR- 1 


ind 

=  ( j -  1 ) *NR  +  i 

r 

=  RI  (i) 

PO 

=  AtmPressure 

P2 

=  S(MvarP,ind) 

P 

=  PO  +  P2 

P 

=  PO 

U 

=  S (MvarU  , ind) 

W 

=  S (MvarW  , ind) 

T 

=  S (MvarT  , ind) 

YK 

=  S (MvarYK, ind) 

Pr 

=  SD (KderPr 

,  ind) 

Pz 

=  SD(KderPz 

,  ind) 

Ur 

=  SD (KderUr 

,  ind) 

Uz 

=  SD (KderUz 

,  ind) 

Urr 

=  SD (KderUrr 

,  ind) 

Uzz 

=  SD ( KderUz z 

,  ind) 

Urz 

=  SD (KderUrz 

,  ind) 

Wr 

=  SD (KderWr 

,  ind) 

Wz  =  SD(KderWz  , ind) 
Wrr  =  SD (KderWrr  , ind) 
Wzz  =  SD (KderWzz  , ind) 
Wrz  =  SD(KderWrz  , ind) 
Tr  =  SD(KderTr  , ind) 
Tz  =  SD ( KderTz  # ind) 
Trr  =  SD(KderTrr  , ind) 
Tzz  =  SD(KderTzz  , ind) 
YKr  =  SD ( KderYKr  ,  ind) 
YKz  =  SD (KderYKz  ,  ind) 
YKrr  =  SD (KderYKrr, ind) 
YKzz  =  SD (KderYKzz , ind) 


Pt 

=  O.dO 

i  set 

to 

zero 

here 

when 

computing 

EQO 

ut 

=  O.dO 

!  set 

to 

zero 

here 

when 

computing 

EQO 

wt 

=  O.dO 

!  set 

to 

zero 

here 

when 

computing 

EQO 

Tt 

=  O.dO 

i  set 

to 

zero 

here 

when 

computing 

EQO 

YKt 

=  O.dO 

i  set 

to 

zero 

here 

when 

computing 

EQO 

WK  =  MolecularWeights 
WM  =  MixtureWeight  (YK,  WK) 
RHO  =  P*WM/ (RU*T) 

PDYN  =  P2 
GZ  =  GRAY 


CP  =  Specif icHeat (T) 
CPt  -  O.dO 

LACP  =  LambdaCp(T) 

LACPt  =  Lambda Cp_ddT(T) 
PRN  =  PrandtlNumber 

LEK  =  LewisNumbers 


Specif icHeat_ddT (T) 
lambda /Cp 
d(lambda/Cp) /dT 


DK  =  LACP/ (LEK*RHO) 

DKt  =  LACPt/ (LEK*RHO) +LACP/ (LEK*RHO*T) 

UC  =  dot_product (DK, YKr )  !  correction  vel.  for  mass  conservation 

UCr  =  dot_product (DK, YKrr) +dot_product ( (DKt*Tr ) , YKr) 

WC  =  dot_product (DK, YKz )  !  correction  vel.  for  mass  conservation 

WCz  =  dot_product (DK, YKzz) +dot_product ( (DKt *Tz ), YKz) 


NUK 
NUKWK 
I  (nu_0  * 
sstoich  = 
Yfuel  = 
Yox  = 

Yfuel_F  = 
Yox_A 
Z 

PHI 
omega 
SOURCE  = 

q 

HEAT 


SignedStoichiometricCoef f s 
NUK*WK 

Wk_0) / (nu_F  *  Wk_F)  @  stoichiometric  conditions 
NUKWK (2) /NUKWK (1) 

YK(1) 

YK  ( 2  ) 

l.dO  !  mass  fraction  of  fuel  in  the  Fuel  Stream 

0.232d0  !  mass  fraction  of  oxydizer  in  the  Air  Stream 
Z_mixfrac ( sstoich,  Yfuel,  Yox,  Yfuel_F,  Yox_A) 
PHI_equivrat io  (Z,  Yfuel_F,  Yox_A) 

MolarProductionRate ( PHI ,  RHO,  YK,  Wk,  T) 

(NUKWK/abs (NUK(l) ) )  *  omega 
HeatReleasePerMole (PHI) 
q  *  omega 


density  gradients,  time  derivative 
RHOrP  =  (Pr*r*RHO*U) /P 


RHOrT  =  - (r*RHO*Tr*U) /T 

RHOrY  =  sum (  ( (r*RHO*U*WM*YKr) /WK)  ) 

RHOr  =  RHOrP+RHOrT+RHOrY 

RHOzP  =  ( Pz*r*RHO*W) /P 

RHOzT  =  - (r*RHO*Tz*W) /T 

RHOzY  =  sum (  ( (r*RHO*W*WM*YKz) /WK)  ) 

RHOz  =  RHOzP+RHOzT+RHOzY 


RHOt 


O.dO  !  leave  off  time  terms  till  later... 


convective  terms 


CONVECT (MeqnP) 
CONVECT (MeqnU) 
CONVECT (MeqnW) 
CONVECT (MeqnT) 
CONVECT (MeqnYK) 


=  (U*RHOr 

=  RHO* (U*Ur 
=  RHO* (U*Wr 
=  RHO* (U*Tr 
=  RHO* (U*YKr 


+ 

W*RHOz ) 

★ 

r 

+ 

W*Uz ) 

★ 

r+*2 

+ 

W*Wz ) 

★ 

r 

+ 

W*Tz ) 

★ 

r 

+ 

W*YKz) 

★ 

r 

if  (iUPWIND. eq. 1  .and.  iDSPACE . eq.  1)  then 

call  UPWIND (ind, i , j , r , RHO, CONVRDER , CONVZDER , CONVECT) \ 
endif 


steady  part  of 

the 

residuals 

EQO (MeqnP  , ind) 

= 

r*RHOt+RHO*U+r*RHO*Ur+r*RHO*Wz+CONVECT (MeqnP) 

EQO (MeqnU  , ind) 

= 

(4  *LACP*PRN*U) / 3 .d0+ (2*LACPt*PRN*r*Tr*U) /3 . dO 

Sc 

& 

- (4*LACP*PRN*r*Ur) / 3 .dO+Pr* (r*r) 

Sc 

& 

- ( 4 *LACPt*PRN*Tr*Ur* (r*r)  )  /3 .dO 

Sc 

& 

- (4*LACP*PRN*Urr* (r*r) ) /3 . dO+RHO*Ut* (r*r) 

Sc 

Sc 

-LACPt*PRN*Tz*Uz* (r*r) -LACP*PRN*Uzz* (r*r) 

Sc 

Sc 

-LACPt*PRN*Tz*Wr* (r*r) 

Sc 

Sc 

- (LACP*PRN*Wrz* (r*r) ) /3 .dO 

Sc 

Sc 

+ (2*LACPt*PRN*Tr*Wz* (r*r) ) /3 .dO 

Sc 

Sc 

+ CONVECT (MeqnU) 

EQO (MeqnW  , ind) 

- 

Pz*r-GZ*r*RHO+ (2*LACPt*PRN*Tz*U) /3 .dO 

Sc 

Sc 

+  ( 2  *LACPt*PRN*r*Tz*Ur ) /3 .dO- (LACP*PRN*r*Urz) /3 . dO 

Sc 

Sc 

- (LACP*PRN*Uz) / 3 . dO-LACPt*PRN*r*Tr*Uz-LACP*PRN*Wr 

Sc 

Sc 

-LACPt*PRN*r*Tr*Wr-LACP*PRN*r*Wrr 

Sc 

Sc 

+  r*RHO*Wt- ( 4  *LACPt *PRN*r*Tz*Wz ) /3 .dO 

Sc 

Sc 

- ( 4  *LACP*PRN*r*Wzz) / 3 . dO 

Sc 

Sc 

+ CONVECT (MeqnW) 

EQO (MeqnT  , ind) 

- 

- ( (HEAT*r ) /CP) -LACP*Tr-LACP*r*Trr+r*RHO*Tt 

Sc 

Sc 

-LACP*r*Tzz 

Sc 

Sc 

- (CPt*LACP*r* (Tr*Tr ) ) /CP-LACPt*r* (Tr*Tr) 

Sc 

Sc 

- (CPt*LACP*r* (Tz*Tz ) ) /CP-LACPt*r* (Tz*Tz) 

Sc 

Sc 

+ CONVECT (MeqnT) 

EQO (MeqnYK, ind) 

- 

- ( r* SOURCE) + RHO*UC * YK+ r * RHOr *UC *YK+r* RHO* UCr*YK 

Sc 

Sc 

+r*RHOz*WC*YK+r*RHO*WCz*YK- (LACP*YKr) /LEK 

Sc 

Sc 

- (LACPt*r*Tr*YKr) /LEK+r*RHO*UC*YKr 

Sc 

&  - (LACP*r*YKrr) /LEK+r*RHO*YKt - (LACPt*r*Tz*YKz) /LEK  & 

&  +r*RHO*WC*YKz- (LACP*r*YKzz ) /LEK  & 

&  +CONVECT (MeqnYK) 

ENDDO 

ENDDO 


RESIDUAL  -  BOUNDARY  NODES 


CONDITIONS  AT  THE  INLET:  BOUNDARY  1 

j  =  1 

DO  i  =  2,  NR 


ind  =  ( j -1) *NR  +  i 


r 

= 

RI  (i) 

PO 

= 

AtmPressure 

P2 

= 

S (MvarP, ind) 

P 

ii 

o 

+ 

fO 

P 

= 

PO 

U 

=  S (MvarU 

,  ind) 

w 

=  S (MvarW 

,  ind) 

T 

=  S (MvarT 

,  ind) 

YK 

=  S (MvarYK 

,  ind) 

Pz 

— 

SD (KderPz 

,  ind) 

Ur 

= 

SD (KderUr 

,  ind) 

Uz 

= 

SD (KderUz 

,  ind) 

Urz 

= 

SD (KderUrz 

,  ind) 

Wr 

= 

SD ( KderWr 

,  ind) 

Wz 

= 

SD (KderWz 

,  ind) 

Wrr 

= 

SD ( KderWrr 

,  ind) 

Wzz 

= 

SD ( KderWz z 

/  ind) 

Tr 

= 

SD (KderTr 

,  ind) 

Tz 

= 

SD ( KderTz 

,  ind) 

Wt  =  O.dO 

if  (iDTIME.gt.l)  & 

Sc  Wt  =  VELinlet (i) *MODULATION*2*PI/PERIOD*cos (timephase) 

WK  =  MolecularWeights 
WM  =  MixtureWeight  <YK,  WK) 

RHO  =  P*WM/ (RU*T) 

PDYN  =  P2 
GZ  =  GRAV 

LACP  =  LambdaCp (T)  !  lambda/Cp 

LACPt  =  LambdaCp_ddT (T)  !  d ( lambda/Cp) /dT 

PRN  =  PrandtlNumber 

Winlet  =  SB1 (MvarW  ,i) 


Tinlet  =  SB1 (MvarT  ,i) 
YKinlet  =  SB1 (MvarYK, i ) 


& 

& 

& 

& 

& 


EQO (MeqnP  , ind) 


EQO  (MeqnU  , ind) 
EQO (MeqnW  , ind) 
EQO (MeqnT  , ind) 
EQO (MeqnYK, ind) 


Pz*r-GZ*r*RHO+ ( 2*LACPt *PRN*Tz*U) /3 .dO  & 

+ (2*LACPt *PRN*r*Tz*Ur) /3 .dO- (LACP*PRN*r*Urz ) /3 . dO  & 
- { LACP*PRN*Uz ) /3 . dO - LACPt *PRN*r*Tr*Uz - LACP*PRN*Wr  & 
-LACPt *PRN*r*Tr*Wr+r*RHO*U*Wr-LACP*PRN*r*Wrr  & 

+r*RHO*Wt- (4*LACPt *PRN*r*Tz*Wz ) /3 . dO+r*RHO*W*Wz  & 
- (4*LACP*PRN*r*Wzz) /3 .dO 
U 

W- Winlet 
T-Tinlet 
YK- YKinlet 


ENDDO 


.Special  treatment 
i  =1 
ind  =  1 
U 

w 

T 

YK 
Pr 

Winlet  = 

Tinlet  = 

YKinlet  = 

EQO (MeqnP 
EQO (MeqnU 
EQO (MeqnW 
EQO (MeqnT 


(for  pressure  BC) 


S (MvarU 
S (MvarW 
S (MvarT 
S (MvarYK 
SD (KderPr 
=  SB1 (MvarW 
=  SB1 (MvarT 


ind) 

ind) 

ind) 

ind) 

,  ind) 

i) 
i) 


SB1 (MvarYK,  i) 


,  ind) 
,  ind) 
,  ind) 
.  ind) 


EQO (MeqnYK, ind) 


Pr 

U 

W-Winlet 
T-Tinlet 
YK- YKinlet 


at  symmetry- inlet  corner  pt 


<<<<<  ...  CUT  ...  >>>>> 


TIME-DERIVATIVES 


iDTIME=0  :  STEADY  STATE  -  NO  TIME  DERIVATIVE  TERMS 


IF  ( iDTIME . eq . 0 )  THEN 

do  ind  =  1, nodes 
do  keq  =  1 , nvar 

F(keq,ind)  =  EQO (keq, ind) 
enddo 
enddo 


goto  999 


iDTIME=l  :  BACKWARD  EULER  DISCRETIZATION  IN  PSEUDO-TIME 


ELSEIF  (iDTIME.eq.l)  THEN 

INTERIOR  POINTS 

do  j  =2 , NZ- 1 
do  i=2,NR-l 


ind  =  ( j-1) *NR+i 

r  =  RI  (i) 

PO  =  AtmPressure 
P2  =  S(MvarP,ind) 

P  =  PO  +  P2 
P  *  PO 

T  =  S (MvarT  , ind) 

YK  =  S (MvarYK  ,  ind) 

WK  -  MolecularWeights 
WM  =  MixtureWeight  ( YK,  WK) 

RHO  =  P*WM/ (RU*T) 

CT ( 2 , ind)  =  RHO  *  r**2 
CT ( 3 , ind)  =  RHO  *  r 
CT ( 4 , ind)  =  RHO  *  r 
do  k  =  1 , NSP 

CT ( 4+k, ind)  =  RHO  *  r 
enddo 

Pt  =  (S (MvarP  , ind) -SI (MvarP  , ind) ) /DTO 

Tt  =  (S (MvarT  , ind) -SI (MvarT  , ind) ) /DTO 

YKt  =  (S (MvarYK, ind) -SI (MvarYK, ind) ) /DTO 

RHOtP  =  ( Pt *r*RHO) /P 

RHOtT  =  - ( r*RHO*Tt ) /T 

RHOtY  =  sum (  ( ( r*RHO*WM*YKt ) /WK)  ) 

RHOt  =  RHOtP  !  +RHOtT+RHOtY 
F (MvarP, ind)  =  RHOt  +  EQO (MvarP, ind) 
do  k  =  2 , NVAR 

dSdt(k,ind)  =  (S (k, ind) -SI (k, ind) ) /DTO 
F(k,ind)  =  CT(k,ind)  *  dSdt(k,ind)  +  EQO (k, ind) 
enddo 

enddo 

enddo 


<<<<<  ...  CUT  ...  >>>>> 


SOME  VARIABLES  MAY  BE  FIXED 


999  continue 

if  (KVFIX.gt.O)  then 
do  ind=l / NODES 
do  k=l , NVAR 

if  (NFIX(k) .eq. 1)  then 

F (k , ind) *S (k, ind) -SX (k, ind) 

EQO (k , ind)  =  0 . dO 
endi  f 
enddo 
enddo 
endif 

!  monitor  norms  of  steady  and  unsteady  residuals 


call  NORM2 (NEL , EQO , eqOnrm) 
call  NORM2(NEL/F/  fnrm) 
dnel=sqrt (dble (NEL) ) 
f nrm=  f nrm/ dnel 
eqOnrm=eqOnrm/dnel 

if  ( ISTEP . ne . -1 )  then 

write (  6,*)  1 @@@  eqOnrm  = 

1 , eqOnrm 

write (  6,*)  ' @@@  fnrm  =  1 

'  ,  fnrm 

write (16,*)  ■@@@  eqOnrm  = 

'  , eqOnrm 

write  (16,*)  ' @@@  fnrm  =  1 

1 ,  fnrm 

endif 

RETURN 

END 

! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
!  %%  %% 

!%%  JACOBIAN  %% 

!  %%  %% 

! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


SUBROUTINE  JACOBIAN 


USE  DimensioningParameters 
USE  ProblemParameters 

USE  MethodParameters ,  ONLY  :  iDSPACE, iDTIME,  & 

Sc  TFLOW ,  DTO  ,  DTI ;  DT2  ,  & 

&  iMATVEC,  Sc 

Sc  KVFIX,  NFIX 

USE  Pointers 

USE  Discret izedProblem  ,  ONLY  :  NVAR ,  NSP,  NR,  NZ,  NODES,  NEL,  & 

Sc  NDRTP ,  Sc 

Sc  RI ,  Z  J ,  Sc 

Sc  SF,  S,  SI,  S2 ,  SX,  Sc 

Sc  SD ,  & 

&  DJ,  DQ,  CTDJ ,  F,  Sc 

Sc  VELinlet 

USE  Differentiation 
USE  OneStepChemistry_Methane 
USE  Simplif iedTransport 


USE  Time_IO_Debug,  ONLY  :  CPUREST,  SECONDM 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

double  precision  YK(NSP_P) ,  YKt (NSP_P) ,  YKinlet (NSP_P) 

double  precision  YKr(NSP_P),  YKz(NSP_P),  YKrr(NSPJP),  YKzz (NSP_P) 

double  precision  LACP,  LACPt ,  LACPtt,  LEK(NSP_P) ,  DK(NSP_P),  DKt (NSP_P) ,  & 
!  stoichiometric  coeffs  (#  mols  produced/destroyed) 

Sc  NUK  (NSP_P)  ,  Sc 

!  mass  coeffs:  NUK(k)*WK(k) 

Sc  NUKWK(NSP_P)  ,  Sc 

Sc  WK(NSP_P),  SOURCE  (NS P_P)  ,  SOURCEp  (NSP_P)  ,  & 

Sc  SOURCEt (NSP_P) ,  SOURCEy  (NSP_P,NSP_P)  ,  HEATy  (NSP_P)  ,  & 

Sc  omegay  (NSP_P) 

double  precision  DSPERT (NODES_P) ,  SHOLD (NVAR_P, NODES_P) 
double  precision  FO (NVAR_P, NODES_P) 

iNumericalDJ  =1  !  1  -  numerical,  0  -  analytical 

timephase  -  2 *PI *TFLOW/Period 


INITIALIZE 


NDJ  =  NVAR_P  *  NVAR_P  *  NOD  E  S_P  !  size  of  array  DJ 
NDQ  =  NVAR_P*LDQ_P*NODES_P  !  size  of  array  DQ 

call  DZERO ( D J , NDJ ) 

call  DZERO (DJN, NDJ)  !  numerical  DJ 


call  DZERO (DQ, NDQ) 
call  DZERO (CTDJ,NEL) 

if  ( iDTIME . eq . 0 )  then 

1  NO  TIME  DERIVATIVES:  steady-state  solve 
elseif  ( iDTIME . eq . 1)  then 

ddt  =  l . DO/DTO  !  pseudo-time  Implicit  Euler 

elseif  (iDTIME . eq . 2 )  then 

ddt=l. DO/DTO  !  time-dependent  Implicit  Euler 

elseif  ( iDTIME . eq . 3 )  then 

ddt=l. DO/DTO  !  time-dependent  Crank-Nicolson 

elseif  ( iDTIME . eq . 4 )  then 

ddt=3 . DO/ (2 . dO*DTO) !  time-dependent  2nd  order  BDF 
else 

write (*,*)  'STOP:  other  time  discrets.  not  yet  implemented' 
stop 
endif 


i _ 

1  NUMERICAL  DJ 

i _ 


IF  ( iNumericalDJ . eq . 1 )  THEN 
call  RESIDUALS  (-1) 


FO  =  F  !  save  initial  residual 

SHOLD  =  S  !  save  current  solution  S 

DO  kvar  =  1,NVAR 


DO  j  -  2 ,  NZ-1 
DO  i  =  2,  NR- 1 

ind  =  ( j - 1 ) *NR+i 

tmp  =  S (kvar, ind) 

DSPERT ( ind)  =  tmp*l.d-8  +  l.d-8 
S (kvar, ind)  =  tmp  +  DSPERT (ind) 
ENDDO 
ENDDO 


call  RESIDUALS  (-1)  !  compute  perturbed  residual 

DO  j  =2,  NZ-1 
DO  i  =  2,  NR- 1 
ind  =  ( j -1 ) *NR+i 

DO  keqn  =  1 , NVAR 

DJ (keqn, kvar, ind) = (F (keqn, ind) -FO (keqn, ind) ) /DSPERT ( ind) 
ENDDO 
ENDDO 
ENDDO 


S  =  SHOLD 


ENDDO 

call  RESIDUALS  (-1) 


ENDIF 


SPATIAL  DERIVATIVES 


if  ( iDSPACE . eq . 1)  then  !  LO  method 
call  SPATI AL_LO  (S) 
else 

call  SPATIAL_CS  (S) 
endif 


JACOBIAN  -  INTERIOR  NODES 


DO  j  *  2 ,  NZ-1 
DO  i  =  2,  NR - 1 


ind  = 

( j -1) *NR  + 

i 

r  = 

RI  (i) 

PO  = 

AtmPressure 

P2  = 

S (MvarP, ind) 

P 

=  PO  +  P2 

P 

PO 

U 

=  S (MvarU  , 

ind) 

w 

=  S (MvarW  , 

ind) 

T 

=  S (MvarT  ; 

ind) 

YK 

=  S  (MvarYK  , 

ind) 

Pr  = 

SD (KderPr 

,  ind) 

Pz  = 

SD (KderPz 

,  ind) 

Ur  = 

SD (KderUr 

,  ind) 

Uz  = 

SD (KderUz 

,  ind) 

Urr  = 

SD (KderUrr 

,  ind) 

Uzz  = 

SD (KderUzz 

,  ind) 

Urz  = 

SD ( KderUr z 

,  ind) 

Wr  = 

SD (KderWr 

,  ind) 

Wz  = 

SD (KderWz 

,  ind) 

Wrr  = 

SD ( KderWr r 

,  ind) 

Wzz  = 

SD ( KderWz z 

,  ind) 

Wrz  = 

SD ( KderWr z 

i  ind) 

Tr  = 

SD ( KderTr 

,  ind) 

Tz  = 

SD ( KderTz 

,  ind) 

Trr  = 

SD ( KderTrr 

,  ind) 

Tzz  = 

SD ( KderTz z 

,  ind) 

YKr  = 

SD (KderYKr 

,  ind) 

YKz  = 

SD ( KderYKz 

,  ind) 

YKrr- 

SD (KderYKrr, ind) 

YKzz  = 

SD( KderYKz z 

,  ind) 

time  ■ 

derivatives 

Pt  =  0 . dO 
Ut  =  O.dO 


YKt  =  O.dO 

IF  (iDTIME.ge.l  .and.  iDTlME . le . 3 )  THEN  !  ps-time  IE,  IE,  CN 

Pt  =  ddt* {S (MvarP  , ind) -SI (MvarP  , ind) ) 

Ut  =  ddt* {S {MvarU  , ind) -SI (MvarU  , ind) ) 

Wt  =  ddt*(S(MvarW  , ind) - SI (MvarW  , ind) ) 

Tt  =  ddt*{S{MvarT  , ind) -SI (MvarT  , ind) ) 

YKt  =  ddt * {S {MvarYK, ind) -SI (MvarYK, ind) ) 

ELSEIF  ( iDTIME . eq . 4 )  THEN  !  BDF2 

Pt  =  0 . 5d0/DT0* { 3 . dO*S {MvarP  , ind) -4 . dO*Sl {MvarP  , ind) 

+1 . dO*S2 {MvarP  , ind) ) 

Ut  =  0 . 5d0/DT0* ( 3 . dO*S (MvarU  , ind) - 4 . dO*Sl {MvarU  , ind) 

+1 . dO*S2 {MvarU  , ind) ) 

Wt  =  0 . 5d0/DT0* { 3 . dO*S (MvarW  , ind) -4 . dO*Sl (MvarW  , ind) 

+1 . dO*S2 (MvarW  , ind) ) 

Tt  =  0 .5d0/DT0* (3 .dO*S (MvarT  , ind) -4 . dO*Sl {MvarT  , ind) 

+1 .dO*S2 (MvarT  , ind) ) 

YKt  =  0 .5d0/DT0* (3 .dO*S (MvarYK, ind) -4 .dO*Sl {MvarYK, ind) 

+1 . dO*S2 (MvarYK  , ind) ) 

ENDIF 


WK  =  MolecularWeights 

WM  =  MixtureWeight  (YK,  WK) 


RHO  = 

P*WM/ (RU*T) 

PDYN  = 

P2 

GZ 

GRAV 

CP 

=  SpecificHeat (T) 

CPt 

=  0  .dO 

CPtt 

=  O.dO 

LACP 

=  LambdaCp(T) 

LACPt 

=  LambdaCp_ddT (T) 

LACPt t 

=  LambdaCp_d2dT2 (T) 

PRN 

=  PrandtlNumber 

LEK 

=  LewisNumbers 

Specif icHeat_ddT (T) 
Specif icHeat_d2dT2 (T) 
lambda /Cp 
d (lambda/Cp) /dT 
d2 (lambda/Cp) /dT2 


DK  =  LACP/ (LEK*RHO) 

DKt  =  LACPt/ (LEK*RHO) +LACP/ (LEK*RHO*T) 

UC  =  dot_product (DK, YKr )  !  correction  vel .  for  mass  conservation 

UCr  =  dot_product (DK, YKrr) +dot_product { (DKt*Tr ) , YKr) 

WC  =  dot_product (DK, YKz)  !  correction  vel.  for  mass  conservation 

WCz  =  dot_product (DK, YKzz) +dot_product ( (DKt*Tz ) , YKz ) 


NUK 

NUKWK 

1  { nu_0 

sstoich 

Yfuel 

Yox 

Yf uel_F 

Yox_A 

Z 

PHI 

omega 

omegap 

omegat 


=  SignedStoichiometricCoef f s 
=  NUK*WK 

*  Wk_0) / (nu_F  *  Wk_F)  @  stoichiometric  conditions 
=  NUKWK (2) /NUKWK (1) 

=  YK { 1 ) 

=  YK ( 2 ) 

=  l.dO  !  mass  fraction  of  fuel  in  the  Fuel  Stream 

=  0.232d0  !  mass  fraction  of  oxydizer  in  the  Air  Stream 
=  Z_mixfrac (sstoich,  Yfuel,  Yox,  Yfuel_F,  Yox_A) 

=  PHI_equivrat io  (Z,  Yfuel_F,  Yox_A) 

=  MolarProductionRate { PHI ,  RHO,  YK,  Wk,  T) 

=  MolarProductionRate_ddP { PHI ,  RHO,  YK,  Wk,  T) 

=  Mol arProduct ionRate_ddT ( PHI ,  RHO,  YK,  Wk,  T) 


omegay  =  MolarProductionRate_ddYK(PHI,  RHO,  YK,  Wk,  T) 

SOURCE  =  (NUKWK/abs (NUK(l) ) )  *  omega 

SOURCEp  =  (NUKWK/abs (NUK(l) ) )  *  omegap 

SOURCEt  =  (NUKWK/abs (NUK(l) ) )  *  omegat 

!  outer  product  of  NUWK  *  (omegay) *T 

SOURCEy  =  spread( (NUKWK/abs (NUK(l) )) ,2,NSP_P  ) *spread (omegay , 1 , NSP  P) 
q  =  HeatReleasePerMole (PHI) 

HEAT  =  q  *  omega 

HEATp  =  q  *  omegap 

HEATt  =  q  *  omegat 

HEATy  =  q  *  omegay 

LOCAL  COMPONENTS  -  MAIN  BLOCK  DIAGONAL 
IF  ( iNumericalDJ . ne . 1)  THEN 
. . . .Analytical  DJ 

DJ (MeqnP  , MvarP  , ind)  =  - ( (r*RHO*Tt) / (P*T) ) + (RHO*U) /P  & 

&  - (r*RHO*Tr*U) / (P*T) + (r*RHO*Ur) /P  & 

&  - (r*RHO*Tz*W) / (P*T) + (r*RHO*Wz) /P  & 

&  +sum (  ( (r*RHO*U*WM*YKr) / (P*WK) )  )  & 

&  +sum (  ( ( r*RHO*WM*YKt ) / (P*WK) )  )  & 

&  +sum(  ( (r*RHO*W*WM*YKz) / (P*WK) )  ) 

DJ (MeqnP  , MvarU  ,  ind)  =  RHO+ (Pr*r*RHO) /P- (r*RHO*Tr) /T  & 

&  +sum(  ( ( r*RHO*WM*YKr ) /WK)  ) 

DJ (MeqnP  ,MvarW  , ind)  =  ( Pz*r*RHO) /P- ( r*RHO*Tz ) /T  & 

&  +sum(  ( (r*RHO*WM*YKz) /WK)  ) 

DJ (MeqnP  , MvarT  , ind)  =  - ( (Pt*r*RHO) / (P*T) ) + (2*r*RHO*Tt) /T**2  & 

&  - (RHO*U) /T- (Pr*r*RHO*U) / (P*T)  & 

&  + (2*r*RH0*Tr*U) /T**2- (r*RHO*Ur) /T  & 

&  - (Pz*r*RHO*W) / (P*T) + (2*r*RH0*Tz*W) /T**2  & 

Sc  -  (r*RHO*Wz) /T  Sc 

Sc  -sum  (  (  (r*RHO*U*WM*YKr)  /  (T*WK)  )  )  & 

Sc  -sum  (  (  (r*RHO*WM*YKt)  /  (T*WK)  )  )  & 

Sc  -sum  (  (  (r*RHO*W*WM*YKz)  /  (T*WK)  )  ) 

DJ (MeqnP  , MvarYK, ind)  =  - ( (Pt*r*RHO*WM) / (P*WK) ) + (r*RHO*Tt*WM) / (T*WK)  & 
Sc  -  (RHO*U*WM)  /WK-  (Pr*r*RHO*U*WM)  /  (P*WK)  & 

Sc  +  (r*RHO*Tr*U*WM)  /  (T*WK)  -  (r*RHO*Ur*WM) /WK  & 

Sc  -  (Pz*r*RHO*W*WM)  /  (P*WK)  Sc 

Sc  +  (r*RHO*Tz*W*WM)  /  (T*WK)  -  (r*RHO*WM*Wz)  /WK  & 

Sc  -  (2*r*RHO*U*YKr*  (WM*WM)  ) /WK**2  & 

&  - (2*r*RHO*YKt* (WM*WM) ) /WK**2  & 

Sc  -  (2*r*RH0*W*YKz*  (WM*WM)  ) /WK**2 

DJ (MeqnU  , MvarP  , ind)  =  (RHO*U*Ur* (r*r) ) /P+ (RHO*Ut* (r*r) ) /P  & 

&  + (RHO*Uz*W* (r*r) ) /P 

DJ (MeqnU  , MvarU  , ind)  =  (4*LACP*PRN) /3 . d0+ (2*LACPt*PRN*r*Tr) /3 . dO  & 

&  +RHO*Ur* (r*r) 

DJ (MeqnU  , MvarW  , ind)  =  RHO*Uz*(r*r) 

DJ (MeqnU  , MvarT  , ind)  =  (4*LACPt*PRN*U) /3 . dO  & 

Sc  +  (2*LACPtt*PRN*r*Tr*U) /3  .  dO  & 

&  - (4*LACPt*PRN*r*Ur) /3 .dO  & 

&  - (4*LACPtt*PRN*Tr*Ur* (r*r) ) /3 .dO  & 

&  - (RHO*U*Ur* (r*r) ) /T  & 

&  - (4*LACPt*PRN*Urr* (r*r) ) /3 .dO  & 

&  - (RHO*Ut* (r*r) ) /T-LACPtt*PRN*Tz*Uz* (r*r)  & 


&  -LACPt*PRN*Uzz* (r*r) - (RHO*Uz*W* (r*r) ) /T  & 

&  -LACPtt*PRN*Tz*Wr* (r*r)  & 

&  - (LACPt*PRN*Wrz* (r*r) ) /3 .do  & 

&  + (2*LACPtt*PRN*Tr*Wz* (r*r) ) /3 .do 

DJ (MeqnU  ,MvarYK,ind)  =  - { (RHO*U*Ur*WM* (r*r) ) /WK)  & 

&  - (RHO*Ut*WM* (r*r) ) /WK- (RHO*Uz*W*WM* (r*r) ) /WK 

DJ (MeqnW  , MvarP  , ind)  =  - ( (GZ*r*RHO) /P) + (r*RHO*U*Wr) /P+ (r*RHO*Wt) /P  & 
&  + (r*RHO*W*Wz) /P 

DJ (MeqnW  ,MvarU  , ind)  =  (2*LACPt*PRN*Tz) /3 ,dO+r*RHO*Wr 

DJ (MeqnW  ,MvarW  , ind)  =  r*RHO*Wz 

DJ (MeqnW  ,MvarT  , ind)  =  (GZ*r*RHO) /T+ (2*LACPtt*PRN*Tz*U) /3 . dO  & 

&  + (2*LACPtt*PRN*r*Tz*Ur) /3 .dO  & 

&  - (LACPt*PRN*r*Urz) /3 .dO- (LACPt*PRN*Uz) /3 .dO  & 

&  -LACPtt*PRN*r*Tr*Uz-LACPt*PRN*Wr  & 

&  -LACPtt*PRN*r*Tr*Wr- (r*RHO*U*Wr) /T  & 

&  -LACPt*PRN*r*Wrr- (r*RHO*Wt) /T  & 

&  - (4*LACPtt*PRN*r*Tz*Wz) /3 .do- (r*RHO*W*Wz) /T  & 

&  - (4 *LACPt*PRN*r*Wzz) /3 . dO 

DJ (MeqnW  ,MvarYK,ind)  =  (GZ*r*RHO*WM) /WK- (r*RHO*U*WM*Wr) /WK  & 

&  - ( r*RHO*WM*Wt ) /WK- (r*RHO*W*WM*Wz) /WK 

DJ (MeqnT  , MvarP  , ind)  =  (r*RHO*Tt) /P+ (r*RHO*Tr*U) /P+ (r*RHO*Tz*W) /P 

DJ (MeqnT  , MvarU  , ind)  =  r*RHO*Tr 

DJ (MeqnT  , MvarW  , ind)  =  r*RHO*Tz 

DJ (MeqnT  ,MvarT  ,  ind)  =  (CPt*HEAT*r ) / CP**2 -LACPt*Tr-LACPt*r*Trr  & 

&  - (r*RHO*Tt) /T-LACPt*r*Tzz- (r*RHO*Tr*U) /T  & 

&  - (r*RHO*Tz*W) /T- (CPtt*LACP*r* (Tr*Tr) ) /CP  & 

&  - (CPt*LACPt*r* (Tr*Tr) ) /CP-LACPtt*r* (Tr*Tr)  & 

&  + (LACP*r* (CPt*CPt) * (Tr*Tr) ) /CP**2  & 

&  - (CPtt*LACP*r* (Tz*Tz) ) /CP  & 

&  - (CPt*LACPt*r* (Tz*Tz) ) /CP-LACPtt*r* (Tz*Tz)  & 

&  + (LACP*r* (CPt*CPt) * (Tz*Tz) ) / CP*  *  2+HEATt 

DJ (MeqnT  ,MvarYK,ind)  =  - ( (r*RHO*Tt*WM) /WK) - (r*RHO*Tr*U*WM) /WK  & 

&  - ( r*RHO*Tz*W*WM) /WK+HEATy 

DJ(MeqnYK, MvarP  , ind)  =  (RHO*UC*YK) /P- (r*RHO*Tr*UC*YK) / (P*T)  & 

&  + (r*RHO*UCr*YK) /P- ( r*RHO*Tz*WC*YK) / (P*T)  & 

&  + (r*RHO*WCz*YK) /P+ (r*RHO*U*YKr) /P  & 

&  + (r*RHO*UC*YKr) /P  & 

&  + ( r*RHO*UC*WM*YK* YKr ) / (P*WK) + (r*RHO*YKt) /P  & 

&  + (r*RHO*W*YKz) /P+ ( r*RHO*WC*YKz ) /P  & 

&  + (r*RHO*WC*WM*YK*YKz) / (P*WK) 

DJ(MeqnYK, MvarU  , ind)  =  r*RHO*YKr 

DJ(MeqnYK, MvarW  , ind)  =  r*RHO*YKz 

DJ (MeqnYK, MvarT  , ind)  =  - ( (RHO*UC*YK) /T) - (Pr*r*RHO*UC*YK) / (P*T)  & 

&  + (2*r*RHO*Tr*UC*YK) /T**2- ( r*RHO*UCr* YK) /T  & 

&  - (Pz*r*RHO*WC*YK) / (P*T)  & 

&  + (2*r*RHO*Tz*WC*YK) /T**2- (r*RHO*WCz*YK) /T  & 

&  - (LACPt*YKr) /LEK- (LACPtt*r*Tr*YKr) /LEK  & 

&  - ( r*RHO*U*YKr ) /T- ( r*RHO*UC*YKr ) /T  & 

&  - (r*RHO*UC*WM*YK*YKr) / (T*WK)  & 

&  - (LACPt*r*YKrr) /LEK- (r*RHO*YKt) /T  & 

&  - (LACPtt*r*Tz*YKz) /LEK- (r*RHO*W*YKz) /T  & 

&  - ( r*RHO*WC*YKz ) /T  & 

&  - (r*RHO*WC*WM*YK*YKz) / (T*WK)  & 

&  - (LACPt*r*YKzz) /LEK 


FORALL  (k=l : NSP) 

!  roughly  correct:  dF(YK)/dYK  is  full  since  F(YK)  contains  RHO,  a  fn  of  all  YK 

DJ (MeqnYK (k) , MvarYK (k) , ind)  =  RHO*UC+ (Pr*r*RHO*UC) /P- (r*RHO*Tr*UC) /T  & 


&  +r*RHO*UCr+ (Pz*r*RHO*WC) /P  & 

&  - (r*RHO*Tz*WC) /T+r*RHO*WCz  & 

&  - (RHO*UC*WM+YK(k) ) /WK(k)  & 

&  - (Pr*r*RHO*UC*WM*YK(k) ) / (P*WK(k) )  & 

&  +  (r*RHO*Tr*UC*WM  +  YK(k) ) / (T*WK(k)  )  & 

&  -  (r*RHO*UCr*WM*YK(k) ) /WK(k)  & 

&  -  (Pz  +  r*RHO*WC*WM  +  YK(k)  ) / (P*WK(k) )  & 

&  +  (r*RHO*Tz*WC*WM*YK(k) ) / (T*WK(k) )  & 

&  -  (r*RHO*WCz*WM*YK(k) ) /WK(k)  & 

&  -  (r*RHO*U*WM  +  YKr (k) ) /WK(k)  & 

&  -  (r*RHO*WM*YKt (k) ) /WK(k)  & 

&  -  (r*RHO*W*WM  +  YKz (k)  ) /WK(k)  & 

&  -  (2*r*RH0*UC*YK (k) *YKr (k)  & 

&  * (WM*WM) ) /WK(k) **2  & 

&  -  (2*r*RHO*WC*YK(k)  *YKz  (k)  & 

&  *  (WM*WM) ) /WK(k) **2 

END  FORALL 

DJ (MeqnYK  , MvarYK  , ind)  =  DJ (MeqnYK  , MvarYK  , ind)  +  SOURCEy 

!  N.B.  SOURCEy  is  a  matrix 


IF  (iDTIME.eq.3)  THEN  !  Crank-Nicolson  only 

DO  k2  =  1,  NVAR 
DO  kl  =  1,  NVAR 

DJ(kl,k2, ind)  =  0 . 5dO*DJ (kl , k2 , ind) 
ENDDO 
ENDDO 


END  IF 

ADD  TIME  DERIVATIVE  TERMS  TO  LOCAL  COMPONENTS  (DIAGONAL  ONLY) 

IF  ( iDTIME . ne . 0 )  THEN  !  omit  time  derivatives  for  steady  problem 

DJ (MeqnP  , MvarP  ,ind)  =  DJ (MeqnP  , MvarP  , ind)  +  (ddt*r*RHO) /P 
IF  ( iDTIME. ne.l)  THEN 

DJ (MeqnP  , MvarT  , ind)  =  DJ (MeqnP  , MvarT  , ind)  - ( (ddt*r*RHO) /T) 

DJ (MeqnP  , MvarYK, ind)  =  DJ (MeqnP  , MvarYK, ind)  +  (ddt*r*RHO*WM) /WK 
END  IF 

DJ (MeqnU  , MvarU  , ind)  =  DJ (MeqnU  , MvarU  , ind)  +  ddt*RHO* (r*r) 

DJ(MeqnW  ,MvarW  , ind)  =  DJ(MeqnW  ,MvarW  , ind)  +  ddt*r*RHO 

DJ (MeqnT  , MvarT  ,ind)  =  DJ(MeqnT  , MvarT  ,ind)  +  ddt*r*RHO 

FORALL  (k=l : NSP) 

DJ(MeqnYK(k) ,MvarYK(k) , ind)  =  DJ (MeqnYK (k) , MvarYK (k) , ind) 

+  ddt*r*RHO 

END  FORALL 


END  IF 


END  Analytical  DJ 
END  IF 


SPATIAL  COMPONENTS 


DQ (MeqnP 

, KderPr 

,  ind) 

DQ (MeqnP 

, KderPz 

,  ind) 

DQ (MeqnP 

, KderUr 

,  ind) 

DQ (MeqnP 

, KderWz 

,  ind) 

DQ (MeqnP 

, KderTr 

,  ind) 

DQ (MeqnP 

, KderTz 

,  ind) 

DQ (MeqnP 

, KderYKr 

,  ind) 

DQ (MeqnP 

, KderYKz 

,  ind) 

DQ (MeqnU 

, KderPr 

,  ind) 

DQ (MeqnU 

, KderUr 

,  ind) 

DQ (MeqnU 

, KderUz 

,  ind) 

DQ (MeqnU 

, KderUrr 

,  ind) 

DQ (MeqnU 

#  KderUz z 

,  ind) 

DQ (MeqnU 

, KderWr 

,  ind) 

DQ (MeqnU 

, KderWz 

,  ind) 

DQ (MeqnU 

, KderWrz 

,  ind) 

DQ (MeqnU 

, KderTr 

,  ind) 

Sc 

Sc 


DQ (MeqnU 

, KderTz 

,  ind) 

DQ (MeqnW 

, KderPz 

,  ind) 

DQ (MeqnW 

, KderUr 

,  ind) 

DQ (MeqnW 

, KderUz 

,  ind) 

DQ (MeqnW 

, KderUrz 

,  ind) 

DQ (MeqnW 

, KderWr 

,  ind) 

DQ (MeqnW 

, KderWz 

,  ind) 

DQ (MeqnW 

, KderWrr 

,  ind) 

DQ (MeqnW 

, KderWz z 

/  ind) 

DQ (MeqnW 

, KderTr 

,  ind) 

DQ (MeqnW 

, KderTz 

,  ind) 

Sc 

Sc 


DQ (MeqnT 

, KderTr 

/  ind) 

DQ (MeqnT 

, KderTz 

,  ind) 

DQ  (MeqnT 

, KderTrr 

,  ind) 

DQ (MeqnT 

, KderTz z 

,  ind) 

DQ (MeqnYK, KderPr 

,  ind) 

DQ (MeqnYK, KderPz 

,  ind) 

DQ (MeqnYK, KderTr 

,  ind) 

DQ (MeqnYK 

, KderTz 

/  ind) 

FORALL  (k=l : NSP ) 
roughly  correct:  dF(YK)/dYK  is 
DQ (MeqnYK(k) , KderYKr (k) 

Sc 

Sc 

DQ (MeqnYK (k) , KderYKz(k) 
Sc 


=  (r*RHO*U)/P 
=  (r*RHO*W)/P 
=  r*RHO 
=  r*RHO 

=  - ( ( r*RHO*U) / T ) 

=  - ( (r*RHO*W) /T) 

=  (r*RHO*U*WM) /WK 
=  (r*RHO*W*WM) /WK 

=  r*r 

=  ( -4*LACP*PRN*r) /3 . dO  & 

- (4  *LACPt*PRN*Tr* (r*r) ) /3 ,dO  +  RHO*U* (r*r) 

=  - (LACPt*PRN*Tz* (r*r) ) +RHO*W* (r*r) 

=  ( -4 *LACP*PRN* (r*r) ) /3 . dO 

=  - (LACP*PRN* (r*r) ) 

=  - (LACPt*PRN*Tz* (r*r) ) 

=  (2*LACPt*PRN*Tr* (r*r) ) / 3 . dO 
=  - (LACP*PRN* (r*r) ) /3 .dO 

=  (2*LACPt*PRN*r*U> /3 .dO  & 

- (4*LACPt*PRN*Ur* (r*r) ) /3 .dO  & 

+ (2*LACPt*PRN*Wz* (r*r) ) /3 . dO 
=  - (LACPt*PRN*Uz* (r*r) ) -LACPt*PRN*Wr* (r*r) 

=  r 

=  (2*LACPt*PRN*r*Tz) /3 .dO 
=  - (LACP*PRN)/3.dO-LACPt*PRN*r*Tr 
=  - (LACP*PRN*r) /3 .dO 

=  - (LACP*PRN) -LACPt*PRN*r*Tr+r*RHO*U 
=  ( -4*LACPt*PRN*r*Tz) /3 . dO+r*RHO*W 

=  - (LACP*PRN*r) 

=  (-4*LACP*PRN*r) /3 .dO 

=  - (LACPt*PRN*r*Uz) -LACPt*PRN*r*Wr 
=  ( 2  *LACPt*PRN*U) /3 . dO  & 

+ (2*LACPt*PRN*r*Ur) /3 .dO  & 

- (4*LACPt*PRN*r*Wz)/3.dO 

=  -LACP- (2*CPt*LACP*r*Tr) /CP-2*LACPt*r*Tr  & 

+r*RHO*U 

=  ( -2  *CPt*LACP*r*Tz ) /CP-2 *LACPt*r*Tz  +  r*RHO*W 

=  - (LACP*r) 

=  - (LACP*r) 

=  (r*RHO*UC*YK) /P 
=  (r*RHO*WC*YK) /P 

=  - ( (r*RHO*UC*YK) /T) - (LACPt*r*YKr) /LEK 
=  - ( ( r*RHO*WC*YK) /T) - (LACPt*r*YKz) /LEK 

full  since  F(YK)  contains  RHO,  a  fn  of  all  YK 
( ind)  =  - (LACP/LEK(k) ) - (LACPt*r*Tr) /LEK(k)  & 

+r*RHO*U+r*RHO*UC  & 

+ ( r*RHO*UC*WM*YK (k) ) /WK(k) 

, ind)  =  - ( (LACPt*r*Tz) /LEK(k) ) +r*RHO*W  & 

+r*RHO*WC+ ( r*RHO*WC*WM*YK (k) ) /WK(k) 

=  - ( (LACP*r) /LEK(k) ) 

=  - ( (LACP*r) /LEK (k) ) 


DQ (MeqnYK (k) , KderYKr r (k) , ind) 
DQ (MeqnYK (k) , KderYKzz (k)  , ind) 
END  FORALL 


IF  (iDTIME . eq. 3)  THEN  !  Crank-Nicolson  only 

DO  k2  =  1,  NVAR*NDRTP 
DO  kl  =  1,  NVAR 

DQ(kl,k2, ind)  =  0. 5dO*DQ (kl , k2 , ind) 
ENDDO 
ENDDO 

ENDIF 

ENDDO 

ENDDO 


JACOBIAN  -  BOUNDARY  NODES 


CONDITIONS  AT  THE  INLET:  BOUNDARY  1 

j  =  1 

DO  i  =  2 ,  NR 


ind 

= 

( j -1) *NR  +  i 

r 

= 

RI  (i) 

P0 

= 

AtmPressure 

P2 

= 

S(MvarP  ,  ind) 

P 

= 

P0  +  P2 

P 

=  P0 

U 

= 

S(MvarU  ,  ind) 

W 

= 

S(MvarW  ,  ind) 

T 

= 

S (MvarT  , ind) 

YK 

= 

S (MvarYK, ind) 

Pr 

= 

SD (KderPr  , ind) 

Pz 

= 

SD(KderPz  , ind) 

Ur 

= 

SD(KderUr  , ind) 

Uz 

= 

SD(KderUz  ,  ind) 

Urr 

= 

SD(KderUrr  , ind) 

Uzz 

SD(KderUzz  # ind) 

Urz 

= 

SD(KderUrz  , ind) 

Wr 

= 

SD(KderWr  , ind) 

Wz 

= 

SD(KderWz  , ind) 

Wrr 

= 

SD ( KderWrr  , ind) 

Wzz 

= 

SD(KderWzz  , ind) 

Wrz 

= 

SD(KderWrz  , ind) 

Tr 

= 

SD (KderTr  , ind) 

Tz 

- 

SD ( KderTz  , ind) 

Trr 

= 

SD ( KderTrr  , ind) 

Tzz 

= 

SD(KderTzz  , ind) 

Sc 


time  derivatives 
Wt  =  O.dO 
if  (iDTIME.gt . 1) 

Sc  Wt  =  VELinlet  (i)  *MODULATION*2*PI/PERIOD*cos  (timephase) 


WK  =  MolecularWeights 

WM  =  MixtureWeight  (YK,  WK) 


RHO  = 

P*MW/ (RU*T) 

PDYN  = 

P2 

GZ 

GRAV 

CP 

=  Specif icHeat (T) 

CPt 

=  O.dO 

CPtt 

=  O.dO 

LACP 

=  LambdaCp(T) 

LACPt 

=  LambdaCp_ddT (T) 

LACPtt 

=  Lambda Cp_d2dT2 (T) 

PRN 

=  PrandtlNumber 

Specif icHeat_ddT (T) 
Specif icHeat_d2dT2 (T) 
lambda/Cp 
d{lambda/Cp) /dT 
d2 (lambda/Cp) /dT2 


DJ  (MeqnP  ,  MvarP  ,  ind)  =  -  (  (GZ*r*RHO) /P)  +  (r*RHO*U*Wr ) /P+ (r*RHO*Wt ) /P  Sc 

Sc  + (r*RHO*W*Wz) /P 

DJ (MeqnP  , MvarU  , ind)  =  (2*LACPt*PRN*Tz) /3 . dO+r*RHO*Wr 
DJ (MeqnP  , MvarW  , ind)  =  r*RHO*Wz 
!  from  time  derivative 

IF  (iDTIME.gt . 1)  Sc 

Sc  DJ  (MeqnP  ,  MvarW  ,  ind)  =  DJ  (MeqnP  ,  MvarW  ,  ind)  Sc 

Sc  +  { r*RHO*ddt ) 

DJ  (MeqnP  ,  MvarT  ,  ind)  =  (GZ*r*RHO)  /T+  (2*LACPtt*PRN*Tz*U)  /3  .  dO  Sc 

Sc  +  (2*LACPtt*PRN*r*Tz*Ur) /3  .dO  Sc 

Sc  -  ( LACPt* PRN* r * Ur z)  /3  .  dO-  ( LACPt * PRN*Uz )  /3  .  dO  Sc 

Sc  - LACPt t*PRN*r*Tr*Uz -LACPt* PRN* Wr  Sc 

Sc  - LACPt t * PRN* r*Tr*Wr-  (r*RHO*U*Wr)  /T  Sc 

Sc  -LACPt *PRN*r*Wrr-  (r*RHO*Wt)  /T  Sc 

Sc  -  (4  *LACPtt*PRN*r*Tz*Wz  )  /3 .dO-  <r*RHO*W*Wz)  /T  Sc 

Sc  -  (4  *LACPt *PRN*r*Wzz )  /3  .  dO 

DJ  (MeqnP  , MvarYK,ind)  =  (GZ* r*RHO*WM)  /WK-  ( r*RHO*U*WM*Wr )  /WK  Sc 

Sc  -  ( r*RHO*WM*Wt )  /WK-  ( r*RHO*W*WM*Wz)  /WK 


DJ (MeqnU  , MvarU  , ind)  =  1 
DJ (MeqnW  , MvarW  , ind)  =  1 
DJ (MeqnT  , MvarT  , ind)  =  1 


FORALL  (k=l : NSP) 

DJ (MeqnYK (k) , MvarYK (k) , ind)  =  1 
END  FORALL 


DQ (MeqnP 
DQ (MeqnP 
DQ (MeqnP 
DQ (MeqnP 
DQ (MeqnP 
DQ (MeqnP 
DQ (MeqnP 
DQ (MeqnP 
DQ (MeqnP 


, KderPz 

,  ind) 

= 

r 

, KderUr 

,  ind) 

= 

( 2 *LACPt *PRN*r*Tz )  /3 . dO 

, KderUz 

,  ind) 

= 

- (LACP* PRN) / 3 . dO-LACPt*PRN*r*Tr 

, KderUrz 

,  ind) 

= 

- (LACP* PRN* r) /3 .dO 

, KderWr 

,  ind) 

= 

- (LACP* PRN) -LACPt*PRN*r*Tr+r*RHO*U 

, KderWz 

,  ind) 

= 

( -4  *LACPt *PRN*r*Tz ) /3 ,dO  +  r*RHO*W 

, KderWrr 

,  ind) 

= 

- ( LACP* PRN* r) 

, KderWzz 

,  ind) 

= 

( -4*LACP*PRN*r) /3 .dO 

, KderTr 

,  ind) 

= 

- ( LACPt* PRN* r*Uz) - LACPt* PRN* r*Wr 

& 

& 


DQ (MeqnP  , KderTz  , ind) 


(2*LACPt*PRN+U) /3 .dO 
+ (2*LACPt*PRN*r*Ur) /3 . dO 
- ( 4  *LACPt *PRN*r*Wz) /3 . dO 


& 

& 


ENDDO 


! . Special  treatment 

at  symmetry- inlet  corner  pt 

i  =1 

ind  =  1 

DJ (MeqnU  , MvarU 

, ind)  =  1 

DJ (MeqnW  ,MvarW 

, ind)  =  1 

DJ (MeqnT  , MvarT 

, ind)  =  1 

FORALL  (k=l : NSP) 

DJ (MeqnYK (k) , MvarYK (k) , ind)  =  1 
END  FORALL 

DQ (MeqnP  , KderPr  ,  ind)  =  1 


<<<<<  ...  CUT  ...  >>>>> 


SOME  VARIABLES  MAY  BE  FIXED 


if  (KVFIX.gt.O)  then 
do  ind=l , NODES 
do  k= 1 , NVAR 

if  (NFIX (k) . eq . 1 )  then 
do  kk=l , NVAR 


kkderind= (kk-1) +NDRTP 

DJ (k , kk , ind)  =  0 . dO 

DQ (k , kkderind+1 , ind)  = 

0 

.dO 

DQ (k , kkderind+2 , ind)  = 

0, 

.dO 

DQ (kr kkderind+3 , ind)  = 

0, 

.dO 

DQ (k, kkderind+4 , ind)  = 

0, 

.dO 

DQ (k, kkderind+5 , ind)  = 

0, 

.dO 

enddo 

DJ (k, k, ind)  =  l.dO 
endif 
enddo 
enddo 
endi  f 


RETURN 

END 


